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CHAPTER 1 


INTRODUCTION 


Among today’s flow-solvers for the Euler and Navier-Stokes equations, many 
are based on upwind differencing. Prominent in use are Godunov-type schemes 
[1], in which the upwind bias is achieved by using the solution to the Riemann 
problem defined at each cell face. Riemann’s initial-vglue problem, which is a 
mathematical representation of the shock-tube problem [2], is well-understood and 
easy to model. A membrane separating a gas at two different states is ruptured, 
and shock, contact, and expansion waves are emitted when the two states interact. 
In a Godunov-type scheme this event is supposed to happen at any cell face and 
any time level. 

The Riemann problem can be solved exactly with an iterative method, as 
Godunov [3] did, or approximately, as Roe [4] did, leading to the concept of 
the “approximate Riemann solver.” In Roe’s method the Euler equations are 
linearized about an average state and then solved exactly. Eigenvectors of the 
averaged flux Jacobian, which represent different types of waves, are allowed to 
propagate with speeds equal to their corresponding eigenvalues. These waves 
describe the difference in states across each cell face. 

The “upwind” direction for each wave is clear in one- dimensional flow: it is 
either forward or backward, according to the sign of each eigenvalue. In two or 
three dimensions the direction of wave propagation is not so straightforward: the 


1 
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waves can travel in infinitely many directions. In most current multidimensional 
upwind flow-solvers, however, the upwinding direction is taken normal to the 
cell faces. Thus the schemes are strongly coupled to the grid on which they are 
implemented. Consequently, high resolution of flowfield discontinuities such as 
shock or shear waves can be achieved only when the discontinuities are aligned 
with grid lines. The Riemann solver interprets such waves incorrectly when they 
lie oblique to the grid; this improper interpretation can lead to smearing in the 
numerical solution. 

In recent years, in an attempt to improve the accuracy of flow solutions, a 
number of researchers have developed multidimensional upwind methods where 
information is obtained from or propagated in directions other than the grid con- 
travariant directions. Initial investigations focused on upwind finite-difference 
schemes for advection-dominated flows, e.g. Raithby [5], Hassan et al. [6], and 
Lillington [7]. In these schemes, attempts were made to convect information in the 
streamwise direction, independent of the grid orientation. Jameson [8,9] developed 
a rotated difference scheme for the transonic potential equation. The equation is 
written in a system of coordinates aligned with and normal to the local streamwise 
direction. When the local flowfield is supersonic, grid-aligned derivatives that are 
used to form derivatives in the streamwise direction are upwind-differenced, while 
all other derivatives are centrally-differenced. 

Moretti [10], and later Verhoff and O’Neil [11], developed nonconservative 
characteristic-based schemes for the Euler equations. These schemes define Rie- 
mann variables, and employ grid-decoupled computational stencils determined 
from the directions of the characteristics. The multidimensional Euler equations 
reduce to a set of ordinary differential equations, coupled only through the source 
terms due to entropy variation, which are typically small. Each ordinary differen- 
tial equation describes the propagation of an individual wave along its characteris- 
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tic. Because the schemes are nonconservative, conservation errors are introduced 
when shock waves are present in the solution. These errors can be eliminated by 
applying a shock correction after each iteration. 

Goorjian [12] extended the method of Jameson for use with the Euler and 
Navier-Stokes equations. This method was later improved by Obayashi and Goor- 
jian [13]. In the latter method, grid-aligned input states are used to solve Riemann 
problems in both the local streamwise and normal directions. In the normal di- 
rection the Riemann problem involves two acoustic waves only. 

Colella [14] designed a predictor-corrector algorithm for systems of hyperbolic 
conservation laws in which the left and right states at each grid face are modified 
by waves traveling parallel to the interface. A standard grid-aligned Riemann 
solver is then used to obtain the flux across the face. This is still a direction-split 
approach, in which an oblique wave would be represented by two grid-aligned 
waves. A similar multidimensional method was developed independently by van 
Leer [15]. 

Davis [16] developed a finite- volume method for the Euler equations in which 
the difference formula and computational stencil vary with angle of assumed wave 
propagation. The angle at each cell face is determined by the velocity- difference 
direction, which is normal to a hypothetical steady shock wave that exists between 
the given states bordering the cell face. Derivatives in the grid-aligned frame 
are written in terms of derivatives in the rotated frame. Fluxes normal to the 
assumed frame are calculated using flux-vector splitting. The flux function along 
the assumed shock is a central-difference flux with an arbitrary parameter that 
insures stability. Davis’ method locates steady shock waves very accurately, but 
is unable to locate steady contact discontinuities. 

Roe [17] designed a multidimensional method based on the decomposition of 
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local gradients. The number and type of waves present is chosen, then angles, 
strengths, and speeds are determined from the local data. One model uses four 
acoustic waves with orientation ninety degrees apart, one entropy wave with ar- 
bitrary orientation, and a shear wave represented by a uniform vorticity. Initial 
efforts to implement this model have been made by Kroner [18] and Struijs et al. 

[19] . 

Hirsch and Lacor [20] decouple the numerical solution from the grid-direction 
by seeking an approximate diagonalization of the Euler equations. The equations 
are written in terms of entropy, a component of velocity, and two acoustic-like 
variables. In general, the Euler equations are not diagonalizable because the two 
Jacobian matrices (in two dimensions) are not simultaneously diagonalizable. In 

[20] , however, the similarity transformation is based on derivatives of the solution, 
allowing more freedom. Still, a source term may arise, which can be minimized. 
Two physical directions play a crucial role in the definition of flow variables and 
characteristic directions. One is aligned with the local pressure gradient, and the 
second is related to the strain-rate tensor. In practice these two directions are 
frozen to improve convergence. A variation of Raithby’s [5] streamline-upwinding 
scheme is used to interpolate the flowfield variables, since standard grid-direction 
interpolation was found to produce no advantage over grid-aligned methods. 

Powell and van Leer [21] and Powell [22] implemented a cell-vertex scheme 
for quadrilateral grids consisting of two basic steps: a residual calculation and a 
residual distribution. In the first step the residual is calculated based on a flux 
integral, and in the second step the residual is sent in a weighted manner only 
to the nodes that define the “downwind” face. The convection directions and 
corresponding characteristic quantities chosen are the same as those derived by 
Hirsch and Lacor [20]. These directions can be frozen to help the convergence of 
the scheme. A similar cell-vertex scheme was instituted by Giles et al. [23]. This 
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method differs from Powell’s primarily in that the computations are performed on 
triangular grids using a trapezoidal rule integration for the residual calculation, 
and the local flow direction (rather than characteristic directions) is used to define 
the “downwind” face. 

Levy et al. [24] developed a method which expands upon the work of Davis 
[16]. A dominant-direction angle is chosen, either in the direction of the local 
pressure gradient or the flow velocity, and two sets of left and right states are 
obtained at each cell face via interpolation from the surrounding flowfield data. 
One of these sets of left and right states is aligned with the dominant-direction, 
while the other set is aligned normal to the first. Then, two Riemann-type solvers 
are used to obtain fluxes in the rotated frame. The components of these fluxes in 
the grid direction are added to obtain the flux at the face. 

Parpia and Michalek [25] independently derived a grid-independent upwind 
finite-volume method for the Euler equations very similar to the method proposed 
in this thesis. Left and right states at an interface are still interpolated along grid 
lines, but a multidimensional four-wave pattern made up of two acoustics, a shear, 
and an entropy wave is assumed to describe the difference in states. The strengths 
of these waves are chosen such that the sum of the jumps in the flow properties 
across the waves is minimized. 

Finally, Dadone and Grossman [26] developed a rotated upwind scheme for 
the Euler equations in which flux-difference splitting is applied along two orthog- 
onal directions for each cell face. The directions are determined from the pressure 
gradient, and the left and right states are selected from appropriately chosen cells 
near to the face where the flux is being computed. 

One of the common features of many of the upwind methods described above 
is that information needed at special points because of physical considerations 
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must be interpolated from surrounding data points. Since this adds an extra level 
of complexity to any method, particularly for non-Cartesian meshes, it was decided 
early in the development of the present scheme to only use information obtained 
by interpolation along grid lines as the input to the approximate Riemann solver. 
It would then be the job of the solver to make “intelligent” use out of information 
gleaned from these left and right states. This constraint puts some limits on the 
ability of the solver to recognize what is going on in the flowfield, but the resulting 
simplicity and low expense of the method seem to outweigh its drawbacks. 

The current method uses five waves to describe the difference in states at a 
grid face. Four of these are acoustic, shear, and entropy waves which act in the 
velocity-difference direction (the same dominant-direction chosen by Davis [16]), 
while the fifth is a shear wave that propagates at a right angle to the other four 
(also used by Parpia and Michalek [25]). This fifth wave allows the method to 
capture oblique steady shear waves sharply. The propagation directions can be 
frozen to improve convergence. The method also makes use of the linearizations of 
the Euler equations due to Roe [4] in order to maintain as simple and inexpensive 
a scheme as possible. 

This thesis is organized as follows. Chapter 2 briefly describes the two- 
dimensional Euler and Navier-Stokes equations in Cartesian and generalized co- 
ordinates, as well as the traveling wave form of the Euler equations. The spatial 
and temporal discretization for both explicit and implicit time-marching schemes 
are described in Chapter 3. Chapter 4 outlines the grid-aligned flux function of 
Roe [4], while Chapter 5 details the derivation of the 5- wave grid-independent 
flux function. Chapters 6 and 7 contain stability and monotonicity analyses of 
the 5-wave model, respectively. Two-dimensional results are provided in Chapter 
8. The extension to three dimensions is made in Chapter 9, with corresponding 
results given in Chapter 10. Finally, Chapter 11 gives conclusions. 



CHAPTER 2 


GOVERNING EQUATIONS 


2.1 Two-Dimensional Navier-Stokes Equations 


The two-dimensional Navier-Stokes equations in Cartesian coordinates can 
be written in conservation form as 


<9U <9F dG 3F V dG v i \ 

H + Bi + 55 = ~dT + 8y' 

where the ~ indicates dimensional variables. The conserved variables are U = 
[p } pu,pv,pE^ T , where p is the density, u and v are the components of velocity in 
the x and y directions, and E is the specific total energy. The inviscid flux vectors 


are 


F = 


pu 

pu 2 + p 

puv 

puH J 


G = 


pv 

puv 

pv 2 + p 
L pvH 


(2.2) 


where the specific total enthalpy H = E +p/p. The viscous fluxes are 


F = 

J. V 


0 

Til 

r 2 i 


VVjTij + fc|? J 


G v — 


0 

Tl2 

T22 


vm + k% j 


(2.3) 


where 


dVi dvA , rdv k 
fii = ^W i + ex i ) + dxJ ii - 


(2.4) 
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In equations (2.3) and ( 2 . 4 ), summation convention is implied and Vi u, V 2 v, 
X\ - x, and X 2 = j/- 

The ideal-gas equation-of-state closes the set of equations: 


P = ( 7 ~ 1)P 




(2.5) 


The above equations can be nondimensionalized as follows. First the Prandtl 
number, Reynolds number, and freestream Mach number are defined by 


Pr = 


PCp Re = hsM M 00 = l^, (2.6) 

Ooo 


k Poo 

where is freestream velocity and l is some characteristic 

length, specific to the problem considered. Then each of the variables is nondi- 
mensionalized via: 


, = *22 , = v i = il M = f- 

£ Poo a oo Mo© 


A Xi P_ 

A = — — X{ — P — - -2 

Poo t PooO-oo 


E 

E = —- 

a* 


(2.7) 

( 2 . 8 ) 


Assuming a calorically perfect gas, T is replaced by the expression 

f = 


a 2 


(7 - l)c p ’ 


(2.9) 


and a is nondimensionalized by a^. 

After substitution, the nondimensional Navier-Stokes equations can be writ- 
ten 


dU dF dG dF v dG 


dt + dx + dy dx " r dy 


+ 


( 2 . 10 ) 


where 


U = 


P 

pu 

pv 

ipE J 


( 2 . 11 ) 
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pu 


’ pv 


F 7 

= 

pu 2 + p 

puv 

G = 

puv 

pv 2 + p 




. puH . 


. pvH . 



r 

0 


0 

- 



Til 

G v = 

T"l2 




T21 


T 22 



VjT\j — Qi . 


-VjT2j -Q 2 j 


( 2 . 12 ) 


(2.13) 


The nondimensional shear-stress and heat-flux terms are given by 

Moo f 




Re rv 

Qi = - 


dVi dVj\ L .«Vfc \ 
dXj + dXi) +X dX k 6x} ) 

(2.14) 

Moo/x d(a 2 ) 

RePr( 7 - 1) dXi ' 

(2.15) 


Again Vx = u, V 2 = v, X x = x, and X 2 = y. The equation-of-state is 

( u 2 + v 2 \ . 

P = ( 7 - l )p y E 2 ) ' 


(2.16) 


The Navier-Stokes equations can also be written in generalized curvilinear 
coordinates, where the coordinate directions are defined by 

( = ({X ’ y) (2.17) 

V = v( x iV)- 

Using the chain rule, the derivatives in curvilinear coordinates are written in terms 
of the derivatives in Cartesian coordinates: 


ri .1 

B( 

_Q_ 

-i 


x i Vi 

[ x v Vv\ 


r JL i 

6x 

_Q_ 

By 


(2.18) 


The determinant of the 2 x 2 matrix in (2.18) is defined as the inverse of the metric 
Jacobian J: 

J~ 1 = x^y,, - x v y£. (2.19) 

Inverting (2.18), one can obtain the following equations for the metric terms: 

(x = Jy v iy = ~ Jx v 
t] x = —Jyt Vy = Jx t- 


( 2 . 20 ) 
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The nondimensional Navier-Stokes equations (2.10) can be written in terms of the 


curvilinear coordinates: 


^ + |({.F + « V G) + |(,,P + % G) = 
^ (£*F V + ^j/G v ) + (t/iF v + T] y G v ) . 


( 2 . 21 ) 


Multiplying (2.21) by J -1 , applying the chain rule, and combining and cancelling 
terms, the two-dimensional Navier-Stokes equations in curvilinear coordinates and 
conservation form become: 


dU* , dF* ( dG 

~dT + ^ + 


dl ?* 


+ 


dGl 


d( dr] di dr, ’ 


where 


U* = - 

J 


F* = — 


pu 

pu*u + £ X P 
pu*v + (yP 
pu* H 


P 

pu 

pv 

L pE 




pv 

pv*u + T) x p 

PV*V + T] y p 

pv* H 


F* = - 

v J 




0 


£»Tli + £y T 12 

(x T 21 + iy T 22 

lUVjTij - Ql) + £y{Vj T 2j - Qi) J 
0 

VxTn + f]yT\2 
Vx T 21 + Vy T 22 
Vx{VjTij - Ql) +Vy( V i T 2 j - Q*)\ 

u* = ixU + iyV 


(2.22) 


(2.23) 


(2.24) 


(2.25) 


(2.26) 


(2.27) 


V* = 7J X U + T] y V. 

The terms u* and v* are the contravariant velocity components, and Vi and V 2 
represent u and v, respectively. T%j and Qi are still given by (2.14) and (2.15), but 


now 


_iL. — £ — + 77 — 

dXi ~^ x dt Vx drj 

d A 

dX 2 ~^di + Vy drj‘ 


(2.28) 
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Therefore, the shear stress and heat flux terms are 


Mqq 
R e 


r u = — - {(A + 2#x)(£,t » i + Vxu v ) + + VyVn)} 


M r 


r 22 = {(A + 2/i)(£,V < + VyVv) + + VxUr,)} 


Re 


Tl2 


= r 21 = ^r(tv u t + Wv + t* v * + W'l) 
«■ = -! 


(2.29) 

(2.30) 

(2.31) 

(2.32) 

(2.33) 


The ideal gas equation of state is still given by (2.16). 

In the Navier-Stokes equations, Stokes’ hypothesis, X + (2/3)/x = 0, is used 
for bulk viscosity. Also, 7 is taken as 1.4 and Pr is taken as 0.72. Sutherland s 
law for molecular viscosity, 

3/2 


- JL - JL 

M A 00 V^o 


Tqq + c 
T + c 


(2.34) 


is e 


mployed, with f* = 460°R and c = Sutherland’s constant = 198.6°R. 


2.2 Traveling- Wave Form of the Euler Equations 

Many numerical methods for the Euler equations, including those discussed 
in Chapters 4 and 5, are based upon the knowledge that certain types of waves 
are emitted when two fluid parcels at different states interact. The directions 
and strengths of these waves define the way in which information is propagated 
through the domain. 

In the upwind-differencing method, care is taken that numerical information is 
propagated similarly. It is therefore instructive to look at traveling-wave solutions 
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to the Euler equations. Upwind-differencing also applies to the convective and 
pressure terms of the Navier-Stokes equations. The viscous terms of the Navier- 
Stokes equations are always centrally-differenced. 

The Euler equations are the same as the Navier-Stokes equations (2.1), (2.10), 
or (2.22), except that the viscous terms F v and G v (or F* and G* ) are taken as 
zero. Starting with equation (2.10), the nondimensional Euler equations in two 
dimensions are written as 


au dF dG_ 

dt dx dy : 


(2.35) 


where U is given by (2.11) and F and G are given by (2.12). These equations can 
also be written in quasilinear form: 


aw aw aw 

— - + A^— + B-r— = 0, 
dt dx dy 


(2.36) 


where W is the vector of primitive variables, W = [p,u,v,p] T , and A and B are 
the matrices 


A = 


B = 


" u 

P 

0 

o - 


0 

u 

0 

i Ip 


0 

0 

u 

0 


.0 

pa‘ 

0 

u 


' V 

0 

P 

o ■ 


0 

V 

0 

0 


0 

0 

V 

i Ip 


.0 

0 

pa 2 

V . 



(2.37) 


(2.38) 


Traveling- wave solutions to (2.36) are of the form 

W {x,y,t) = W(xcos0 + ysinfl - At), 


(2.39) 


where 6 is the angle that defines the direction of wave propagation. Insertion of 
(2.39) into (2.36) results in the eigenvalue problem 


(Acos0 + Bsin#)t>W^ — AtfW, 


(2.40) 
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where $W is the amplitude of the traveling wave. The four eigenvalues and 

corresponding right eigenvectors yielded by (2.40) are: 

A 4 = ucosQ + vsin # + a 


A 2 = ucosd + usin# — a 
A 3 = ucosd + vsin# 


(2.41) 


A 4 = ucosd + usin# 


Pi 

P 2 

P 3 


a 


a 


1, — cos#, — sin#, a 

P P 


1, — -cos#, — -sin d,a z 

P P 


a 


0, sin#, — cos#, 0 

P P 


(2.42) 


P 4 = [1,0,0, 0] T . 

These eigenvectors represent: (1) an acoustic disturbance that propagates with 
speed \i, (2) an acoustic disturbance that propagates with speed A 2 , (3) a shear 
wave, and (4) an entropy wave. The latter two waves travel with speed A 3 = A 4 , 
the projection of the fluid velocity in the direction of wave propagation. Together 
the four eigenvectors form the matrix P . 

The characteristic variables for the quasilinear form of the Euler equations 
can be computed from 


6W* = P -1 5W, 


(2.43) 


where 


p-i 


’0 

0 

0 

1 


pcos$ 

2a 

pcosd 

2a 

_ p sing 

a 

o 


gain# 

2 a 

_ PsinO 
2a 
pcosfl 
a 

o 



This gives 


6 W* 


l{££ + (cos#f#u + sin#f#v)}‘ 

f {|f - (cos#£#u + sin#£#u)} . 

— ^sin##u + £ cos dSv 
° 

&P ~ I? 


(2.44) 


(2.45) 
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Rewriting the Euler equations in terms of the characteristic variables yields 


dW* 

dt 


+ P -1 AP 


aw* 

dx 


+ P 1 BP 


aw* 

dy 


= 0 . 


(2.46) 


In general, A and B do not commute, so a single 9 that simultaneously diag- 
onalizes both matrices cannot be found. This is indicative of the fact that in 
multidimensional flow the waves propagate in infinitely many directions. As nu- 
merical schemes are limited to modeling the flow with a finite number of waves, the 
choice of wave type and direction of propagation is not trivial. The more “physi- 
cally relevant” the wave types and directions of propagation are, the more likely 
the model will be able to resolve a wide variety of flow features accurately. As 
will be discussed in Chapter 4, grid-aligned wave models choose the grid-normal 
direction as the direction of wave propagation. The grid-independent model de- 
rived in Chapter 5 allows waves to travel in directions dictated by the physics of 

the local flowfield. 


CHAPTER 3 


SPATIAL AND TEMPORAL DISCRETIZATION 


3.1 Finite Volume Formulation 


The Navier-Stokes equations are cast in finite-volume form, which is a dis- 
cretization of the integral form of the equations. Equation (2.22) is integrated 
over a computational cell of area A: 

J j fEl d A + JJ {(f* - F;) e + (g* - g;m = o. (3.1) 

The second integral is converted to a line integral over the boundary 5 of the cell 
using Gauss’ Theorem, 


JJ ™l dA+ J {(F* - F* v )dr, - (G* - G ■;)<£} = 0 . 


(3.2) 


The first integral can be interpreted as the rate of change of the vector of averaged 
conserved quantities (2.23) within each cell. The factor 1/ J can be taken outside 
of the time derivative when the computational grid is fixed with respect to time. 
The line integral is discretized over the four faces of each cell; in the present 
formulation it is assumed that each cell is a quadrilateral. The final result is the 
finite- volume form of the Navier-Stokes equations. 


= -(S^F* + 6 v G*)i,j + {6tK + s n G v)i,i 

J dt 

= s»,j> 


(3.3) 
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where the right-hand side residual term S is made up of an inviscid and a viscous 
part: S — Sj “j - Sy ^ and 

{*<(•)}« = — . 

(3.4) 

The terms (F and (G*)ij±i/2 represent the inviscid fluxes normal to the 
cell faces. These fluxes are evaluated at the cell faces through the use of a flux 
function , which is the primary concern of the present study. A grid-aligned flux 
function is discussed in Chapter 4, and a grid-independent flux function is derived 
and discussed in Chapter 5. The second parenthetical expression in (3.3) consists 
of viscous fluxes at the cell faces, which are determined via central differencing. 

The following relationships exist between the Jacobian J , the metric deriva- 
tives, and the cell areas and cell face lengths in the finite-volume formulation: 


1/J = Cell area A 

^j/(JAs^) = i component (i = x or y) of unit normal to £=constant cell face 
of length 

T]i/( JAsjj) = i component (z = x or y) of unit normal to 77 =constant cell face 
of length A 3 V . 


Using these relations, it is sometimes convenient to write (3.3) in a different form: 


dU,,j 

dt 


~r~ ( E ~ 

Ai 'j l/=i i=i 




(3.5) 




P9g 

# _ P9g u + pcos0 g 
pq g v + psin0 9 
pq g H 

0 

Tucosfly + T 12 sin0 ff 
t 21 cos# 9 + T 22 sin0 9 

(Vj^j - Qi)cos0 9 + (VjT 2 j - Q 2 ) sin0 9 _ 


(3.6) 


(3.7) 
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where Q g is the angle that the outward-pointing cell face normal makes with the 
x-axis, and q g is the outward velocity normal to the cell face, given by 

q g = ucos 9 g + vsin 6 g . (3-8) 

Furthermore, $*A s t is the inviscid normal flux at cell face t, evaluated through 
the use of a flux function. The term ($ V )/As/ is the viscous normal flux at cell 
face evaluated using central differencing. 


3.2 Explicit Time-Marching 


The solution for the vector of averaged conserved variables can be advanced 
explicitly in time using an m-stage time-stepping scheme. Using the definition of 
S from (3.3), the scheme is: 

U (1) = U (n) + 771 JAfS (n) 

D (2) = u( n) +T7 2 JAiS {1) 

: (3.9) 


U(m-i) = D (n) +7?m _ 1 j A ts< m - 2) 

UCn+i) = Xj(n) + JAtS (m_l) , 

where the superscript n denotes the current time level and n + 1 denotes the next 
time level. The superscripts l,2,...,m - 1 denote intermediate time levels or 
stages. The coefficients 1 7l , 772 , • - Vm - 1 are chosen to give desirable damping 
and stability properties for the scheme. When m = 1 the scheme reduces to the 
one-step forward-Euler time-stepping scheme 


U("+!) = U (n) + JA/S (n) . 


(3.10) 
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3.3 Implicit Time-Marching 


The solution can also be advanced implicitly in time. Starting with (3.3), the 
left-hand side is discretized to first-order in time and the right-hand side terms 
are taken at time level n + 1 : 


1 A t U< n ) 
J At 


-(« { f* + h, G*)<" +1 > + (»<f; + s, 


(3.11) 


where A t U (n) = U (n+1) - U (n) . The right-hand side terms are linearized about 


time level n: 



(3.12) 


(3.13) 


Notice that in (3.13) the viscous matrix Jacobian terms are split in two parts: a 
matrix with derivatives that are a function of { only and a matrix with derivatives 
a function of 77 only. 

The spatial cross- derivative terms ^[F*]u(* 7 ) and ^[G*]u(0 are treated 
explicitly, lagged in time, while ^[F*]u(0 and 6 v [G*) v {v) as well as the inviscid 
matrix Jacobian terms are treated implicitly. Equation (3.11) becomes: 


I 

JAt 


+ 61 


/of* 
V dV 


sf; 

dU 




= - (X(F* - F*) + ^ n (G* - g;)) 


(3.14) 


=T. 
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The equations are approximately factored and solved for A t U^ n ^ in two sweeps 


I f (d F 
J At + 8( \ dU 

dG* dG 


+ fir, 


j At 77 V au au 


) 


A t U' = T 


( 3 . 15 ) 


A ‘ U(n) =JA^ A ‘ U ' 


where the term A t U' is an intermediate result. The conserved variables then are 
updated at the cell centers using 


U(n+i) _ u(n) + A t U (n) . 


( 3 . 16 ) 


The implicit spatial derivatives of the convective and pressure terms are spatially 
first-order accurate, resulting in block tridiagonal inversions for each sweep. For 
example, the left-hand side of the first sweep in (3.15) is a block banded matrix 
with the following structure for the ith row: 


. . . , 0, — A+ 1 , (A+ * - A- 4 + 1/( JAf )) , A- 4 , o, 


( 3 . 17 ) 


where represents the portion of (aF*/aU — dF*/5U(£)) at ce ^ face * + 2 

contributed from the left (the zth cell), and A7 +1/2 represents the portion con- 
tributed from the right (the (z + l)st cell). Each of these terms can be divided 
into inviscid and viscous parts: 

A+ = (aY +(A v r 

»+, V )i+x \ /»+i ( 318 ) 

A >~+i = ( Ai ) i+i + ( Av ) i+ j- 

Since the viscous terms are centrally differenced, the viscous Jacobians are given 


by 



( 3 . 19 ) 
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However, using the flux functions discussed in Chapters 4 and 5, it is very diffi- 
cult and computationally tedious to obtain exact expressions for the inviscid flux 
Jacobians. Instead the following approximations are used: 


/ \+ lf/0F*\ \dF* \ 

[ Ai )i + i “2\UuJ i + |9U j+ J 

( Y -if (—\ — \ 

t A, h + }“2\lauA +1 su i+ J 


(3.20) 


The second terms within the braces are Roe-averaged terms at the interfaces, 


equal to 


= R A R"\ 


where R and A. are defined in Chapter 4, along with the Roe-averaged (hatted) 


variables. 

The left-hand side of the second sweep of (3.15) is of similar form to (3.17). 
However, the A’s are replaced by B’s, which are functions of derivatives of G* 
and G* . The approximate inviscid Jacobians for the second sweep are given by 



dG* 1 

av Y 

e g- \ 

“ ™ i+i ) 


(3.22) 


The appropriateness of using (3.20) and (3.22) for the left-hand side inviscid Ja- 
cobians when the grid-independent model is employed on the right-hand side is 
discussed within the context of stability in Chapter 6. 



CHAPTER 4 


GRID-ALIGNED FLUX FUNCTION 


Most flux computations have two distinct stages: a projection stage and an 
evolution stage. In the projection stage of a finite-volume scheme, left and right 
states are obtained at interfaces via interpolation along grid lines from surrounding 
cell-center values. For first-order spatial differencing, the state variables (normally 
the primitive variables) are extrapolated to a cell face k + ^ (where k represents 
the grid index i or j of a structured 2-D grid) using 


W L = W fc 


W R = w* +1 . 


(4.1) 


For higher-order spatial differencing, 

W L = W fc + i[(l - k) A_ + (1 + k)A+]* 

Wr = W fc+1 - i[(l - k)A+ + (1 + «)A_] fc+1 , 


(4.2) 


where 


(A + )j. = W*j +1 — W*. 


(4.3) 


(A_)* = Wfc — Wfc-!. 

When k = — 1 ? (4.2) gives second-order fully upwind spatial differencing, while 
k = 1/3 gives third-order upwind-biased differencing. Limiting of higher-order 
terms can be employed at this stage of the grid-aligned flux function in order to 
eliminate numerically-induced oscillations near regions of high gradient such as 
shock waves. 
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In the evolution stage of the flux computation, the flux at the interface is 
computed as a function of the left and right states obtained from (4.1) or (4.2). 
As discussed in the introduction, in a Godunov-type solver this stage numerically 
models the physical process defined by the one-dimensional Riemann problem. 
This process is illustrated in figure 4.1. At time zero, a membrane separating the 
left and right states ruptures, and a shock wave, a contact discontinuity, and an 
expansion fan propagate into either side, with strengths and velocities depending 
upon the initial conditions. The flux at the interface can be determined when 
these strengths and velocities are known. 

The grid-aligned solver of Roe [4] is an approximate Riemann solver, in which 
the Euler equations are linearized about an average state and solved exactly. When 
used in a two-dimensional scheme the eigenvectors of the matrix in the linearized 
system of equations, representing acoustic, shear, and entropy waves, are assumed 
to propagate in a direction normal to the grid interface. 

The grid-aligned flux function of Roe, representing the inviscid flux * (3.6) in 
the two-dimensional Euler equations or Navier-Stokes equations, is given below. 
The flux at a face is computed using any one of the following three equations (all 

three are equivalent): 


$ = $i, + 'y ^ Aj.f2j.Rfc 

A<0 

$ = ^ AfcflfcRfc 

A>0 

1 4 

<& = - (#l + ^r.) _ ~ y^iAfc|f2feRfe- 

2 z k= 1 


(4.4) 


The last equation can be interpreted as a central difference term plus a dissipation 
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(4.5) 


term. The eigenvectors are given by 

* 1 T 

Ri = [l, u + acos 8g,v + asin9 g ,H + aq g \ 

R 2 = [l,u — acos9 g ,v — asm9 g ,H — aq g j 

R-3 = [0, — asin0 ff , dcos^j, df ff ] 

R 4 = [l,u,v, i(d 2 +d 2 )] T 

for the equations written in conserved-variable form. These eigenvectors corre- 
spond to the eigenvectors Pfc given by (2.42) for the equations written in primitive- 
variable form, and represent, respectively, -(-acoustic, —acoustic, shear, and en- 
tropy waves. The kth. wave of this system has a strength fijb, evaluated as the fcth 
component of the vector f 2: 


n = 


r 2p (Ap + pdAg 9 ) 

^(Ap-paAqg) 
\p Ar ff 

L ~gs (d 2 A p — A p) 


(4.6) 


where A(-) = (-)h — (-)l and 


The wavespeeds are 


q g — ucos9 g + vsin9 g 
r g = — usin 9 g +vcos9 g . 

Ai = q g a 
A 2 = q g Oj 
A 3 = q g 
A4 = 9g* 

The Roe-averaged values (denoted by hats) are defined as 

P = sfPLpR 

U = UlW + Uft(l — w) 

V — VLW + Vh(1 — w ) 

H = H l w + H r { 1-w) 


(4.7) 


(4.8) 


(4.9) 


«= \ Kit ~ !) 


H - i (d 2 +t) 2 ) 
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where w = y/pZ/{y/pL + yfpR)- The values q g and f g in (4.5) and (4.8) are given 
by (4.7) with Roe-averaged values u and v replacing u and v. 

The wavestrengths f Ik given in (4.6) result from satisfying the equation 

4 

AU = U R - U L = ^^k- ( 4 - 10 ) 

k=l 

In other words, the sum of the eigenvectors times their corresponding strengths 
describes the difference in states across an interface. A geometric interpretation 
of this will be given below. Since the R* are eigenvectors of the matrix d$/dU, 
the additional equality 

4 

A# = #R - = ^2^k&k&k (4T1) 

k= 1 

is also satisfied. It can be seen from (4.4) that if all of the wavespeeds are positive 
in the grid-normal direction, the flux computed at the interface will be the flux 
from the left, $ L . Conversely, if all the wavespeeds are negative, the flux will be 
computed as $ R . In both cases this amounts to the upwind choice for #. 

The grid-aligned model can be interpreted in a geometric sense by looking at 
the effects of the acoustic and shear waves in (Ati, Av, Ap)-space. (The entropy 
wave only causes a change in the density, so it is not representable in this space.) 
Grid-aligned acoustic waves cause a change in velocity in the grid-normal direction, 
along with a proportional change in the pressure and density according to the 


relations: 


+acoustic : 
—acoustic : 


Sp 

pa 2 

Sp 

pa 2 


Su Sv Sp 

acosOg asin0 g p 

Su Sv _ Sp 

acos0 g asinffg p 


(4.12) 


These expresions can be derived easily from the acoustic wave eigenvectors P x and 
P 2 in (2.42) with 8 taken as 8 g . Figure 4.2(a) depicts the assumed propagation 
direction of the + and - acoustic waves at an arbitrary grid interface, the normal 
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Figure 4.1: The Riemann Problem 



Figure 4.2: Grid-Aligned Acoustic Waves 
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of which makes an angle 6 g with the x-axis, while figure 4.2(b) shows the effect of 
each of the waves in ( Au, Ad, Ap)-space. The change in density A p is not pictured. 
In this second figure, the state at a given point in space changes by an amount 
(Au, Ad, A p) as drawn by the heavy solid lines when a + or - acoustic wave 
passes. The length of the lines in state space are representative of the strengths 

of the waves. 

The shear wave also propagates in the grid-normal direction, as shown in 
figure 4.3(a). However, shear waves cause a change in velocity normal to the 
direction in which they propagate, with no change in the pressure or density; this 
is dictated by the relation: 


shear : 


6u 

as'mOg 


6 v 

acos&g ' 


(4.13) 


This expression can be derived torn P, in (2.42). Hence Ihe velocity change across 
a shear wave is in the direction shown in figure 4.3(b) in (An, An, Ap)-space, where 
again the length of the line in this space is proportional to the strength of the 


wave. 

Given a left state L and a right state R at an interface, the grid-aligned 
flux function of Roe interprets the difference with a combination of +acoustic, 
-acoustic, shear, and entropy waves such that (4.10) is satisfied. An example 
is drawn in figure 4.4. L, the representation of the state to the left of the cell 
face, is placed at the origin, and the right state R is located at (Au,Ad,A p), 
as determined by the differences between L and R. All waves propagate in the 
e g -direction, represented by the vertical plane in the figure. The effects of the two 
acoustic waves and the shear wave are shown in the figure as heavy solid lines (the 
entropy wave is not represented). The acoustic waves cause a change in velocity 
in the ^-direction along with proportional changes in pressure, while the shear 
wave causes a change in velocity normal to B g . Since this is a linearized model, 



y 



b) Effect in (Au, Av, Ap)- Space 


Figure 4.3: Grid- Aligned Shear Wave 



Figure 4.4: Grid-Aligned Wave Decomposition 
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the order in which the waves are taken is immaterial, and there is no difference 
between shock-type and expansion-type acoustic waves. 

A steady shock wave that is aligned with a mesh interface is interpreted 
correctly by this model, as depicted in figure 4.5(a): the difference between the 
states L and R is described essentially by a single acoustic wave. However, the 
grid-aligned method smears a shock wave that lies oblique to the mesh. The 
difference in states in this case cannot be described by a single acoustic wave 
since the velocity-difference vector F fl - V L is not in the ^-direction. Hence the 
model must introduce both a shear wave as well as an acoustic wave of the opposite 
family to account for the discrepancy in AF. These extra waves, depicted in figure 
4.5(b), add dissipation which smears the numerical solution. The cone delineated 
by the dashed lines in figure 4.5(b) is defined by the effects of all acoustic waves 
of a given strength and arbitrary orientation, with one endstate at L. It is referred 
to in Chapter 5 as the “acoustic cone.” 

The grid-aligned flux function can also misinterpret a pure shear wave that 
lies oblique to a grid face. This situation is illustrated in figures 4.6(a) and (b). 
In figure 4.6(a) left and right states are indicated on a (Au, Av, Ap)-diagram. 
There is no pressure difference between L and R, and the velocity-difference vector 
V R - Vl is at some angle other than 90° to the ^-direction (it would be normal 
to the -direction for the case of a shear wave aligned with the gnd face, and 
the wave model would intepret the difference with a single shear wave of the type 
shown in figure 4.3(b)). The grid-aligned scheme now includes two acoustic waves 
in its interpretation of the difference in states. These waves add dissipation which 
smears the numerical solution. Additionally, if the wavespeeds associated with 
each of these acoustic waves are of opposite sign, then the scheme computes a 
flux at the interface with a pressure that is different from the correct pressure by 
an amount A p, as shown in figure 4.6(b). In this figure, a time history of the 
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b) Oblique to Grid 


Figure 4.5: Grid- Aligned Wave Model Interpretation of a Single Shock 
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b) Wave Directions for Subsonic Flow 

Figure 4.6: Grid- Aligned Wave Model Interpretation of an Oblique Shear 

Wave 
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wave locations is drawn in relation to the grid-normal direction. The flux at the 
interface is computed as either the left flux plus the change across left-running 
waves or, equivalently, the right flux minus the change across the right-running 
waves. In either case it can be seen that the incorrect pressure is given at the 
interface. 



CHAPTER 5 


GRID-INDEPENDENT FLUX FUNCTION 

The motivation behind the development of the present grid-independent ap- 
proximate Riemann solver is the desire to be able to recognize and appropriately 
model both shock and shear waves regardless of their orientation with respect to 
the grid. The method for accomplishing this goal is described in this chapter. 

The projection stage of the flux computation is identical to that of the grid- 
aligned method described in Chapter 4. In other words, primitive variables are 
interpolated along grid lines using either (4.1) or (4.2). This stage is different 
from that of many of the grid-independent methods under development by other 
researchers, which assign values to faces via interpolation in some grid-independent 
direction. However, since this latter type of interpolation can be very complicated 
and costly, it was decided early in the development of the current scheme only to 
use grid-aligned interpolation. 

It is then the job of the flux function to make use of this information in an 
intelligent fashion during the evolution stage of the flux computation. This is 
accomplished in the following way: 

(1) A primary direction of wave propagation is chosen that is more repre- 
sentative of the physics of the local flowfield than the direction defined by the 
grid. 

(2) The difference in L and R states is represented with a combination of 
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acoustic, shear, and entropy waves. 

(3) A flux is formed in the grid-normal direction from the information prop- 
agating in the physically-relevant grid-independent direction. 

The next three sections of this chapter will describe in detail the methods 
chosen to satisfy each of these three aspects. 


5.1 Wave Propagation Direction 


The primary wave propagation direction used at each interface is the velocity- 
difference direction 

° d = tan “ 1 (ti) ’ (5,1) 

defined from -f to £. This represents the angle that the velocity-difference vector, 
AV, makes with the x-axis, as shown in figure 5.1 for two arbitrary states V L and 
V R . The ^-direction is chosen because in this frame the velocity components vi 
and vr normal to 0 d are equal, as depicted in figure 5.2. Therefore the differences 
between the two states can be interpreted either as a compression normal to 0 d or a 
shear aligned with B d . In figure 5.3(a), the former interpretation is illustrated. The 
velocity components tangential to the shock are equal (only the normal component 
is affected by the shock). Also, the shock wave could be propagating with some 
velocity u s in the ^-direction. The value of u s is zero for a steady shock wave. 
A shear-wave interpretation of the difference in velocities is illustrated in figure 
5.3(b). Here, the shear wave propagates with velocity v L = vr in the ( 0 d + f )* 
direction. This propagation velocity is zero for a steady shear wave. 

Other choices for the dominant wave- propagation direction, as used by other 
researchers ( e.g . [24]) for grid-independent models include the pressure-gradient 
direction and the flow direction. The first of these is not used in the present 
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investigation because its calculation requires data from surrounding cell centers, 
and not just the left and right states interpolated along grid lines. The flow 
direction is not used because it is not normal to a shock that lies oblique to the 
flow. Hence the resolution is not as crisp as with the use of the velocity-difference 
direction. 


5.2 Wave Decomposition 


Since there are two interpretations (figures 5.3(a) and (b)) of a velocity dif- 
ference in terms of a dominant wave, it necessary that the method be able to 
model both types of waves as well as have some way of determining which is a 
better description of the true situation. The present method models both types 
by describing the difference in states by a combination of two acoustic waves and 
an entropy wave propagating in the ^-direction, and an additional shear wave 
propagating in the (0 d + f )-direction. This shear wave causes a change in velocity 
parallel to 9 d with no change in pressure, thus allowing for sharp capturing of 
oblique shear waves of the type depicted in figure 5.3(b). The propagation and 
effect of this (9 d + f ) shear wave is shown in figures 5.4(a) and (b). The propa- 
gation and effect of the + and -9 d acoustic waves is shown in figures 5.5(a) and 
(b) for comparison. 

The representation in primitive- variable form of the (9 d + f ) shear wave is 


P 3 in (2.42) with 9 = 9 d + 


P (0 d + ir /2)sh.ear 


0 


0 

— pSin (0 d + f) 


— -COS0<1 

p cos (0 d + J ) 


— -sin#<f 

0 


. P 0 


(5.2) 


The two acoustic waves and the entropy wave are represented by Pi, P 2 ) and P 4 
in (2.42), with 9 = 9 d . 




a) Direction of Propagation b ) Effect in (Au, Av, Ap)-Space 

Figure 5.4: (6 d + f ) Shear Wave 



Figure 5.5*. Acoustic Waves 
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The strengths of the four waves must satisfy 

AW = £M\, ( 5 - 3 ) 

Jfe=l 


or, equivalently, 


AU = 53 f 1 *®-*’ 


(5.4) 


Jb=l 


where the R fc represent the waves for the conserved- variable form of the equations: 

- ^T 


Ri = [l,u + acosOd , v + dsin^j, H + aq d } 

A rp 

R 2 = [l,u - acos 6 d , v - asinOd, H - aq d ] 

T 

R 3 = [O,-acos0d, -asin^d, -aqd] 

R,= [1,4, i ,i(i 2 +» ! )] T - 


(5.5) 


The hatted variables are still Roe-averaged variables defined by (4.9), and q d = 
ucos Q d + usin 6 d . Unlike in the grid-aligned method described in Chapter 4, there 
is not a unique combination of these four waves that satisfies (5.4). Although the 
entropy wave always has a strength of 


n 4 = ^r(a 2 Ap- Ap), 

a 1 


(5.6) 


(the same as in (4.6) for the grid-aligned model), there is some freedom in 
picking the strengths of the other three waves. This reflects, as mentioned earlier, 
that there are two types of dominant waves, represented by figures 5.3(a) and 
5.3(b), that could describe the difference in states. The model must choose which 
type is more likely to be representative of the true situation, and allow that type 
of wave to dominate in the numerical representation. 

Two methods that allow the model to choose the “correct” wave type are 
described here. Both are a function of the pressure difference across the interface: 
if a large pressure difference exists, it is more likely that an acoustic wave is 
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primarily responsible for the difference in states. Similarly, a small difference in 
pressure indicates that a shear wave more likely is the primary wave. 

The first method is termed the minimum-pathlength model, and is imple- 
mented by choosing the combination of waves such that the pathlength in (Au, 
Au, Ap)-space is minimized. This minimum-pathlength model is accomplished by 
using either two acoustic waves and an entropy wave or one acoustic, a (9 d + j ) 
shear, and an entropy wave. (Recall that the entropy wave is not representable 
in (Au, Av, Ap)-space.) The choice depends on the location in phase space of the 
right state R relative to the cone defined by all acoustic waves emanating from 
L. By definition, R lies in the 0 d -plane. If R resides inside the “acoustic cone,” 
as is the case with Rj in figure 5.6, then two acoustic waves describe the shortest 
path. If R resides outside the cone, as represented by R 2 in the figure, then one 
acoustic and a (Oj + f) shear wave describe the shortest path. The mathematical 
conditions for R inside or outside the acoustic cone are: 

Inside : (Ap) 2 ^ [pa(Aucos0d + Ausinfl^)] (5*7) 

Outside : (Ap) 2 < [/5d(Aucos5^ + Avsin^d)] 2 . (5*8) 

The minimum-pathlength model always uses three waves out of a choice of four 
possible wave types to describe the difference in states. 

A second strategy is to choose the strengths of the acoustic and shear waves 
such that the path is in some sense closest to the straight line connecting L and R 
in phase space. More specifically, the area between the waves (taken in a certain 
order) and the direct path L-R is minimized. This minimum-area model is due 
to Parpia [27]. A geometric representation is given in figure 5.7, where again the 
entropy wave, although present, is not pictured. If R lies inside the acoustic cone, 
like Ri in the figure, then the path that minimizes the area (shaded region) is 
accomplished by two acoustic waves. If R lies outside the cone, as represented by 
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Figure 5.6: Minimum- Pathlength Model Wave Decomposition (Entropy 

Wave Not Pictured) 


4 



Figure 5.7: Minimum-Area Model Wave Decomposition (Entropy Wave 

Not Pictured) 
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R 2 , then some combination of two acoustic waves and a (9 d + f ) shear wave gives 
the minimum area. The exact expression will be given in the next section and is 
derived in Appendix A. The minimum-area model uses either three or four waves 
to describe the difference in states at each interface. 

Numerical experiments show that both these models can produce nonlinear 
feedback that results in oscillatory flowfields. Small changes in the computed 
values of 9 d feed back into the solution, producing further changes in 6 d . An 
example of this is given in figures 5.8(a) and (b). This is an Euler computation 
of supersonic flow over an airfoil, to be discussed in greater detail in Chapter 8. 
Figure 5.8(a) shows pressure contours over the airfoil using the grid-independent 
model 500 iterations after a restart from a converged solution with the gnd-aligned 
scheme. Velocity vectors showing the ^-directions on the £-faces for the last 
iteration are given in figure 5.8(b). These directions are very oscillatory, causing 
the unrealistic behavior seen in the first figure. In spite of the wild oscillations in 
the flowfield, this solution converges. 

An easy way to inhibit this feedback is to freeze the computed values of 6 d 
at each face at some point in the computation, calling these 9 d . The four wave 
vectors in (5.5) remain the same, only with 0' d replacing 6 d . Since the state R 
does not necessarily he in the *'-plane, at least one additional wave is needed to 
describe the differences between the left and right states at each face. A shear 
wave propagating in the 0 ' d -direction produces a change in velocity normal to 0' d , 
and can therefore be used as the additional wave. It is represented by 

R 5 = [0, — asin^, acos0 d , a(—usin0 d + vcos0' d )] T . (5-9) 

A sketch of this wave is included along with the other four waves using the 
minimum-area model in a (Au, Au, Ap)-space diagram in figure 5.9. When the 
projection of the right state onto the 0' -plane lies within the acoustic cone, then 



a) Pressure Contours 


Figure 5.8: 



b) Od -Directions 


Inviscid Computation of Supersonic Flow Over an Airfoil 
With Od Computed Each Iteration 




Figure 5.9: 5 - Wave Model Wave Decomposition (Entropy Wave Not Pic 
tured 


a) Pressure Contours 



b) ^-Directions 

Figure 5.10: Inviscid Computation of Supersonic Flow Over an Airfoil 

With 6' d Frozen 
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the waves are chosen as shown in the figure for right state R, with projection onto 
the 0^-plane R[. The second case, for which the projection of the right state lies 
outside the acoustic cone, is shown for right state R 2 with projection onto the 
0^-plane R' 2 ; the wave decomposition in the 0^-plane follows the minimum-area 

rule. 

The same airfoil computation is shown in figures 5.10(a) and (b), where the 
e d angles are now computed during the first iteration of restart from a converged 
grid-aligned solution, then frozen for the remainder of the computation. Results 
are now relatively free from oscillations. 

It should be noted at this point that a first-order interpolation procedure is 
used to obtain the left and right velocity values employed in equation (5.1) for 0 d 
even when the overall computation is second-order accurate. This is done since 
first-order interpolation yields smoother variations in 0 d throughout the compu- 
tational domain, giving generally better solution quality. As a consequence, in a 
second-order computation the left and right states L and R obtained via second- 
order interpolation do not necessarily lie in the 0 d -plane, and the fifth shear wave 
R 5 described here is necessary even when the wave propagation directions are not 

frozen. 

5.3 Flux Formulation 

5.3.1 Standard Formula 


The combination of the four waves from (5.5) (with 0 d taken as 9' d ) plus 
the e' d shear wave (5.9) results in a 5-wave model, which generates a family of 
flux formulas with a free parameter 0. This family includes both the minimum- 
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pathlength and minimum-area models discussed above. The flux per unit face- 
length normal to each grid face is calculated using 

* = \ (*l + *r) - \ (5.10) 

2 1 fc=l 

where the five waves are given by 

A p 

Ri = [l,u + acos9' d ,v + asm9' d ,H + aq' d ] 

R 2 = [l,u — acos9' d ,v - asin9' d ,H - aq' d ] 

R 3 = [0, — acos#j, — asint^, —aq' d ] T ( 5 -H) 

R 4 = [l,u,v,^(u 2 +v 2 )] T 
R 5 = [0, — asm9' d ,acos9' d ,ar d ] . 

These represent, respectively: +9' d acoustic, — 9 d acoustic, (9' d + f) shear, 9 d 
entropy, and 9 d shear waves. Also, 


q d = ucos 9 d + wsin 9 d 
r' d = —us\n9 d + vcos9 d . 

The wavestrengths are defined as 

r + p A ( Aucos ^ + Avs{nd d) ■■ 

^ - /3^ (Aucos9' d + Avsin9' d ) 

O = | (/3 — 1)4 (Aucostf^ + Avsin 9 d ) 
a? (a 2 Ap - Ap) 

| (— Ausin 9' d + Aucos 9 d ) 

The minimum-pathlength model is obtained when (3 is taken as 


( 5 . 12 ) 


( 5 . 13 ) 


f3 — min 


Ap/(pa) 

Aucos9' d + AusintfJj ’ 


( 5 . 14 ) 


In this case, the 5-wave model uses only four waves at a time since either the 
(Q' d -|_ shear wavestrength or one of the acoustic wavestrengths is identically 
zero. The minimum-area model results when 


(3 - min 


f A p/(pa) 

\ Aucos 9' d + Avsin#^ 



( 5 . 15 ) 
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and the 5-wave model uses all five waves when R lies outside the acoustic cone, 
and four waves when R lies inside the acoustic cone (the (9' d + f ) shear strength 
is zero). These two expressions are derived in Appendix A. In practice, a small 
number e(=lx 10 -5 ) is added to the denominators in (5.14) and (5.15) to avoid 
division by zero in regions of null gradient. Also, /? is generally limited to be no 
less than 0.05, and is frozen along with 0' d as an aid to convergence. 

Numerical experiments indicate that both the minimum-pathlength and the 
minimum-area models give very similar results, although the minimum-area model 
tends to be slightly more dissipative for a wider range of test cases. Hence, 
it exhibits less oscillatory behavior and usually converges slightly faster. The 
minimum-area model is used for all the computations in Chapters 8 and 10. 

The wavespeed associated with each of the waves in (5.11) is the component 
of the average flowspeed in the direction of wave propagation, plus or minus the 
average speed of sound for the acoustic waves. Since the flux (5.10) is in the 
grid-normal direction, however, it is necessary to take the components of these 
wavespeeds in the 0 g -direction. They are. 

Ax = (<ld + a)cos(6 d — 0 g ) 

^2 = {<id n)cos(fl J j — 0 g ) 

A, = - <y} ( 5 - 16 ) 

A 4 = foot'd - 0g) 

A 5 = Q'd cos (^d — ^ 9 )* 

Notice that this 5-wave model reduces to the grid-aligned approximate Riemann 
solver when 0' d = 6 g and /? = 1 (t.e. the ( 6' d + f ) shear-strength vanishes). 
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5-3.2 Rotated Riemann Solver 

Since the 5-wave model’s wave vectors R*. = Rfc(^d) are not the eigenvectors 
of the matrix d&/dV unless 8' d = 6 g and 0 = 1, it is generally not possible to 

satisfy equation (4.11), i.e., 

5 

Jt=i 

Some question naturally arises, then, as to the appropriateness of equation (5.10) 
for determining the flux. A general form for the flux function is 


$ = 6 

d - 'y 'j AfcflfcRfc 

+ (!-«) 



A fc <o 


>*> 

3 r 

V 

0 

1 


with (0 < 6 < 1). This expression reduces to formula (5.10) for 6 = 0.5, but 
gives different results for # using different 8 since the expressions in brackets are 
not equal in value. This is in contrast to the grid-aligned wave model described 
in Chapter 4, for which (5.18) would give identical * for all 8, as a result of the 
validity of equation (4.11). 

For example, if all A’s at a particular grid face are positive and 8 = 1, the 
flux formula yields $ L , the flux of the state on the “left” of the face. In the 
grid-aligned method this would be the correct answer. If the waves are traveling 
obliquely to the grid direction, however, this is not the answer produced by the 
standard formula (5.10). It may be argued that equation (5.10) is the correct 
formula because of its symmetry - it favors neither input state since it corresponds 
w ith 8 = 0.5 in (5.18) - but this is not always entirely clear. Hence it is difficult 
to say what value of 6, if any, is appropriate in (5.18) for the grid-independent 

model. , 

In order to gain further insight, it is instructive to look at a slightly different 
formulation of the model. In this method, as before, the first step is to determine 
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9 d via equation (5.1), and the second step is to use the five waves (5.11) to describe 
the difference in states. The third step, the determination of the flux normal to 
the cell face, is performed in a different manner, however. The strengths and 
orientations of the five waves are used to define four new states, as shown in figure 
5.11, which lie to the left and right of the face in the 6' d and ( 9 d + -|)-directions. 
Two Roe-type approximate Riemann solvers are then used, one in each direction 
(9 d and 9 d + j), and the resulting fluxes are combined to give a flux on the grid 
face via: 

$ = ^^cos (9 d - 9 g ) - #^ + ^/ 2 sin(^ - 9 g ) (5.19) 

This method is hereinafter referred to as the state-determined model. 

The state-determined model is similar to that of Levy [24] and Dadone and 
Grossman [26]. However, they use information from the surrounding flowfield to 
obtain the four states; the present method uses only the left and right states L 
and R, and applies the information in the acoustic, shear, and entropy waves to 
determine local gradients. 

The four states U L ', Ur', U l ", and Ur" of figure 5.11 are determined as 
follows. First, assume that cos {9-9 g ) > 0 and sin (6-0 g ) > 0. Results are similar 
for other cases. Second, assume for simplicity that 9 d is not frozen, i.e. 9 d = 9 d , 
and that the difference in states is described by the minimum-pathlength model. 
Although the left and right states are interpolated to be at the same location on 
either side of a cell face, an arbitrary finite distance is assumed to separate them 
as shown in the figure. The states Ul' and Ur' are found on the wave fronts 
through L and R of the 9 d waves, while Ul" and Ur" lie on the fronts through 
L and R of the (9 d + f ) wave. These fronts form a rectangle; the states are 
assumed to be in the midpoints of the rectangle’s sides. Furthermore, the waves 
that describe the difference in states are assumed to bridge this difference in a 
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linear fashion. Thus, the 6 d waves have half of their fuU effect on U L " and Ur , 
and the ( 6 d + f ) wave has half of its effect on U L ' and Ur'. 

If equation (5.7) is true, then the difference in states is described by two 
acoustic waves and an entropy wave, all acting in the ^-direction. Hence Ul 
Ul, Ur = Ur, and 

u L " = ur" = u L + \ Yj (5,20) 

fc=l,2,4 

Alternatively, if equation (5.8) is true, then one acoustic and an entropy wave 
in the 0,,-direction and a shear wave in the (6 d + f )-direction are used. In this 


case, the formulas are: 

Ul' = Ul + ^HsRs ( 5 - 21 ) 

Ur' = Ul + nior2Rlor2 + 64R4 + ^sR-3 ( 5 ‘ 22 ) 

Ul" = Ul + ^ (fiio^Ri or 2 + 64R4) + ^Rs ( 5 - 23 ) 

Ur" = Ul + \ (rilor2R.lor2 + 64R4) • ( 5 ‘ 24 ) 


The subscripts (1 or 2) indicate that the acoustic wave R.i or R 2 is used, depending 
upon which minimizes the pathlength in (Au, At,, A P )-s P ace. Roe-averaged values 
are determined once at each face from states L and R, then used in each of the 
Riemann-solvers. 

The state-determined model is naturally fairly expensive, requiring about six 
times as much CPU time as the grid-aligned method to compute the fluxes at the 
faces. It can be shown that the flux based on this state-determined model is not 
the same as the flux (5.18), regardless of the value of 6. However, it has been 
determined numerically that the difference between the two methods is always a 
minimum in the mass flux when 6 = 0.5. Also, if the jump between the left and 
right states is relatively small, then the difference between the two methods in the 
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Figure 5.11: Graphical Representation of State-Determined Model 



a) Pressure-Dominated Difference b) Velocity-Dominated Difference 

Figure 5.12: Typical Plots of Difference Between Standard Flux Formula 

and State-Determined Model 
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other three equations is also usually a minimum near 6 = 0.5. For example, figure 
5.12(a) shows the absolute value of the difference between the fluxes determined 
by the two methods as a function of 6, where the left-state primitive variables 
(p,u,v,p) are (1,2, 3, 4), and the right-state variables are (2,3,4,10). 9 g is taken 
as 0°. These particular states have a relatively large pressure difference, so two 
acoustics and an entropy wave are used to describe the difference in states. A 
second example is given for left and right states of (1,1, 1,1) and (1,2, 1,1) with 6 g 
taken as 45° in figure 5.12(b). In this case, the difference in states is only in the 
velocity, and one acoustic, one shear, and one entropy wave are used to model the 
difference. 

Hence, equation (5.18) with 6 = 0.5 is, in a sense, “best” for this particular 
form of the flux function in that it generally yields a flux very close to that given by 
the state-determined model. In fact, in practice the flowfields produced by the two 
methods are often virtually identical. Results in the remainder of the paper are 
therefore obtained solely with the standard formula for the 5-wave model (5.10) 
(which is the same as (5.18) with 6 = 0.5), since its expense is significantly lower 
than that of the state-determined model. 


5.4 Recovering the Grid-Aligned Scheme 


One unresolved issue facing the 5-wave model is the fact that the model does 
not reduce to the standard grid-aligned model when 8' d = 9 g if (3 1. In other 

words, if the ( 6' d + f ) shear wave has any strength other than zero, the standard 
grid-aligned flux is not obtained when the waves are assumed to travel in the 
grid-normal direction. Intuitively, it seems that this flux should be obtained. For 
example, if all wave speeds are positive and act in the ^-direction, then it seems 



54 


desirable that the grid-aligned upwind flux $l should be computed as the flux 
at the interface. But unless (3 = 1 (along with 9' d = 9 g ), the inequality in (5.17) 
holds and the flux function (5.10) does not yield $ L as a result. 

Also, it is difficult to imagine a steady-state circumstance when 6 d - 9 g and 
a (^ + f ) shear wave exists between the states. As an example, two pictures 
are drawn in figure 5.13 of neighboring cells with a velocity-difference direction 
that coincides with the grid-normal direction. In figure 5.13(a) a relatively large 
pressure difference exists between the states, and a 0' d acoustic wave is presumed 
by the model to describe the difference. In figure 5.13(b), no pressure difference 
exists, and the model describes the difference with a (*i + f ) shear wave. However, 
this shear wave is aligned such that it passes exactly through each cell center. 
Even if these are assumed to lie infinitely close to either side of the wave, this is 
an unlikely situation. This interpretation also is obviously inconsistent with the 
notion that the data in each cell represent cell averages , for with this shear-wave 
orientation the averages in the left and right cell should be identical. 

Unfortunately, it is not clear how to recover the grid-aligned method auto- 
matically when 9' d = 9 g , while still retaining the capability of resolving oblique 
shear waves. Two different methods which have been attempted to resolve this 
issue are described in this section. The first method involves making (3 a function 
of the difference 9' d -9 g . When 9' d = 9 g , the grid-aligned scheme must be recov- 
ered, so (3 must be made equal to 1.0. This can be enforced by choosing a new 0* 

as 

(3* = (3 + 1-/3) {cos [2 (9 d - 0 9 )] + 1} , ( 5 - 25 ) 

where 0 is computed using (5.14) or (5.15). The resulting 0\ used in place of 0 
in (5.13), varies smoothly between 0 and 1.0, depending on the difference between 

0 d and 9 g . 
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There are two significant problems with this method. First of all, resolution 
of pure oblique shear waves is reduced since 0* is near zero only within a narrow 
range of angles ( 8' d - 6 g ). Hence most pure shear waves are now interpreted using 
some acoustic waves as well, resulting in added dissipation and smearing of the 
numerical solution. The second problem arises in certain circumstances, and most 
notably in the case of the computation of flow over an airfoil. Around the airfoil s 
leading edge the grid-normal angles of a structured-grid vary rapidly through 
about 180°. Since the 0^-directions do not vary as much (see figure 5.10(b)), 
the 0* values resulting from (5.25) generally vary rapidly between 0 and 1.0. 
This oscillatory behavior of 0* can cause unrealistic results with “kinks” and/or 
oscillations in contours of flowfield variables near the leading edge. 

Example solutions using 0 and 0* are given in figures 5.14(a) and (b). Shown 
are pressure contours resulting from an Euler computation over a NACA 0012 air- 
foil at M = 0.3, a = 1°, on a 65 x 25 C-mesh. Figure 5.14(a) shows contours using 
0 (without equation (5.25)). Results are fairly smooth, except for irregularities 
due to grid coarseness. Results using 0*, shown in figure 5.14(b), shows two re- 
gions of local maximum pressure at the nose and large “kinks” in the contours in 
front of the nose. One effect of this irregular behavior is the nonconvergence of 
lift to the correct value as the grid is refined. 

Several variations in equation (5.25) have also been attempted, including 
the use of functions which very rapidly transition from 0* = 1 to 0* = 0 near 
\e' d - 9g\ = 0 and ± 7 T, and remain flat at 0 over much of the 8 d - 6 g range. This 
allows for sharper resolution of most oblique shear waves, but the second problem 
relating to airfoil flow remains. 

The second attempted method for recovering the grid-aligned scheme when 
0' d = 8 g involves replacing the (8' d -\-^) shear wave with -j-0 d acoustic, — 8'J acoustic, 
and 6 d shear waves. The 0'J angle is a function of 8 d and 8 g : it is set equal to 


b) Using /r 

Figure 5.14: Pressure Contours Near Leading Edge of NACA 0012 Airfoil 

M = 0.3, a = 1°, 5- Wave Model 
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6' d when 0' d = 9 g , and is smoothly transitioned to 0'J = 0 d + f otherwise. Hence, 
when 6' d = 9 g the grid-aligned scheme is recovered, and when 0 d = 9' d + y the 
5-wave model is recovered. One way to accomplish the transition is to use the 
relation 

0'J = 0' d + imn (2 \0' d - 0 g | , |) . (5.26) 

The strengths of the three new waves which replace the ( 8 d + f ) shear wave are 

= ~ — — t(Aucos 9 d + Avsin0^)(cos0jcos<?y + sin^sin#^) 

d 2 CL 

Q-0'Ja.c = ( 5 . 27 ) 

d « 

£l 0 » shear = (/? — 1 )t(Aucos0j + Ausin0j)(cos0jsin0j — sin0jcos0j), 

* a 

The corresponding components of the wavespeeds in the grid-normal direction are 
^+0"ac = (ucos 0'J + vsin0 d + a)cos(8 d - 8 g ) 

A_ e »ac = (ucos0 d + vsinflj — a)cos(0j - 0 g ) (5.28) 

A#»shear = (uCOS0 d + Vsinfl^COS^ - 9 g ). 

Aside from the obvious drawback that this model requires a total of seven waves 
rather than five, it also (just as the method based on /?*) gives very unrealistic 
results near the leading edge of airfoils. By far the best solutions in general are 
still obtained using the original 5- wave model, in spite of the fact that it does not 
recover the grid-aligned method. 


CHAPTER 6 


STABILITY ANALYSIS 

The stability analyses of both explicit and implicit time-marching schemes 
with the 5-wave model for the Euler equations are discussed in this chapter. Be- 
cause of the complexity of the Navier-Stokes equations, it is much more difficult 
to obtain expressions for its stability. However, it has been found empirically that 
the stability condition for the Navier-Stokes equations for flows with relatively 
high cell Reynolds numbers is usually only slightly more restrictive than the CFL 
condition for the Euler equations [28]. 


6.1 Explicit Time-Marching 


Writing the equations in the finite-volume form (3.5) with viscous terms ig- 
nored, the stability analysis requires the eigenvalues of the Fourier transform of 
the right-hand side of 


SU At A 

At— = — — 


Ot 


/= 1 


( 6 . 1 ) 


= JAtSi. 


The right-hand side can be expanded using the i and j indices corresponding to 
the £ and rj directions of a structured grid: 


Ai ir = + (* Aj )i./+* - (* Ai, )i.i-i) • < 6 - 2 > 
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The normal flux $ per unit face-length at each of the faces is obtained with the 
flux function described in Chapter 5: 


*‘+i = 5{ (#L)t+ i + (#r) ‘+= ) “ (6 3) 

= ^{($l)* + i +(^R) fe+ i} - \ |D| fc+i {(U R ) fc+ i -(U L ) fc+ i|, 

where k represents the index i or j . The following operator symbols are introduced: 


Aj.+ a(’) = (-)*+i “ (-)fc 
MO = (0*+i _ (')fc-i- 


(6.4) 


When equation (6.3) is inserted into (6.2) at the four faces, terms such as A fc+1 / 2 # 
and can arise (depending on the spatial order of accuracy). These terms are 
linearized using 







8 k $ 



(6.5) 


and all d$/dU and |D| terms are also linearized about the cell center using Roe- 
averaged values (e.g. | jt-t-i /2 = l^lfc— 1/2 = |D|). ^ * s further assumed that all 

cells are square with face length As. Linearized, equation (6.2) becomes 


At 


dV 

at 


= -LU, 


( 6 . 6 ) 


where, for first-order spatial differencing 


L 





+KHk’))sj 


(6.T) 
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For second-order fully one-sided spatial differencing. 



The 0 9 *^ and 0 9 P are the directions normal to the grid faces in the i and j directions, 
respectively. In the present analysis, 0 9 * } is prescribed, and 0 9 j) is taken as 0 { g ] + 
tt/2. 

For completeness, the matrix d4/dU is given here: d$/dU = 


0 

( c z)g < fi~'uq g 

(c y ) g 4>-vq g 
{ 2cf)-jE)q g 


( c *)a 

(cr) 9 (2-7 )u+q g 
(c x ) g v-(c y ) g (-f-l)u 

C 42 


( c »)a 

( c s)a** — ( c * (”T 1)^ 
( c »)s( 2 -7)v+g tf 
C 43 


0 

M,{ 7-1) 

(c») 9 (7-1) ’ 


79a 


(6.9) 


where C A2 = (c x ) g (~(E-<t>) - ( 7 -lH and C 43 = {c y ) g {~iE-<t>) - ( 7 -!)^- Also » 


q g = (c x ) fl ii + (c„) 9 t>. 


( 6 . 10 ) 


The terms ( c x ) g and (c v ) 9 represent the components of the grid-face unit normal 
in the x and y coordinate directions, respectively. In other words, ( c x ) g = cos0 9 
and (cy )p — sin0 9 . 

|D| is the matrix that satisfies 


D 



( 6 . 11 ) 


where the summation is over all the waves. For the 5-wave model, the elements 
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of |D| are given by: 


i i 


A 4 


(R 4 )j + 4>i 1 - izu - iiv 


j2 


J3 


= -(7 - i)w 6 + 6 

-(7 - 1)»6 + 6 
= (7 - 1)61 


(6.12) 


where 


{|Ai| (Ri); + A 2 (R 2 )j — 2 |A 4 j (R.4)j| 

£ 2 = ^cos 9' d {|Aj (Ri)j - A 2 (Ra),-} + ~ t cos6' d A 3 (R-s)i 


2a 2 

: £ t 

2d 

- rsin^ A 5 (Rs)j 

a 


(6.13) 


=^rin«i { |Ai | (ft-Oi - \>*\ (Aj)j} + 


(*>)> 


+ tcos^ A 5 (R 5 ); 


and (Ri)j represents the jth element of the vector Ri, as defined in (5.11). The A’s 
are defined in (5.16). It should be noted that for the grid-aligned model described 
in Chapter 4, the |D| matrix which satisfies (6.11) can be written as 

R -1 , (6.14) 


D 



= R 

A 


(grid— aligned) 

dU 




whereas this is not true for the 5-wave model, since its R* are not the eigenvectors 
of the d&/dU matrix. 

In order to obtain the Fourier transform of the right-hand side of (6.6), the 
Fourier symbols for the difference operators are inserted into the equation. These 
symbols are: 


3 


5 


9(A t+i ) = e‘<“’ - 1 


a 


(S) = 


- .«<*> _ „-«'*> 


(6.15) 
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where k represents the index i or j, and and are the wavelengths of 
the perturbations in the i and j directions, respectively. For first-order spatial 
differencing (L defined by (6.7)), the Fourier transform of the right-hand side of 
(6.6) is: 



i(||(^ 0 )) si nC (0 + |d(^>)| (1 - cos^>)+ ( 6 - 16 ) 

where the variable v is the CFL number, defined as 

^={u,(*W)+u,(flO))}|^, (6.17) 

for square cells with face length As. Here, o»(flf ) ) and w(^ j) ) are the maximum 
wave speeds \q\ + a in each of the grid directions. 5(-L) is a complex-valued 4x4 
matrix. 

For stability, the locus of the eigenvalues of the Fourier transform, often called 
the “Fourier footprint,” must lie inside the stability boundary of the time- marching 
scheme. In general, the Fourier footprint of the 5- wave model is a function of fc', 
the Mach number M = y/u 2 + v 2 j a, the flow angle a = tan 1 (v/u), /3, &* d , and dg \ 
as well as the perturbation wavelengths and Since the Fourier footprint 

is a function of so many variables, it is difficult to perform a thorough numerical 
analysis. However, an extensive number of variations in the independent variables 
have been tested. In each case and are both cycled through 17 values from 
— 7 r to 7 r inclusive, and four eigenvalues are obtained at each of the 289 conditions. 
Based on the results obtained the following trends are noted: 

(1) The relative magnitudes of the eigenvalues are strongly dependent on the 
Mach number. In general, increasing M increases the magnitudes, but at larger 
and larger M an asymptotic limit is reached. 
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(2) Many of the modes of the Fourier footprint can have eigenvalues that lie 
on the imaginary axis. In particular, this occurs when a and 6 d differ by 90°. 

From this analysis, it is clear that the explicit forward-Euler time stepping 
scheme would be unstable for many of the modes, since the stability boundary for 
forward-Euler is a unit circle centered at (—1,0) and does not include any part of 
the imaginary axis except the origin. However, 2-stage or higher schemes can be 
designed that satisfy this requirement. For example, the 4-stage scheme 

u (1) = U (n) + T/JAtSi (n) 

Tj( 2 ) = U (n) + - JAtS| (1) 

? (6.18) 

Xj( 3 ) = U (n) + - JAtSi (2) 

2 

u (n+1) = U (n) + JAtSi (3) • 

has a stability region including a finite part of the imaginary axis whenever 77 < 
0.6756. Its amplification factor is given by P(z) = 1 + z + ^z 2 + jz 3 + qV 2 * • The 
stability boundary is defined by |P(z)| — 1- 

An attempt was made to devise a “worst case” Fourier footprint for this 
scheme by choosing independent variables that yield the largest eigenvalue extent 
in the Real- Imaginary plane. Of all the combinations of variables tested, the one 
that gave the largest footprint was: M = 100, a = 90°, (3 = 0, 8g^ = 0°, and 
0' d = 22.5°. Given this footprint, an 77 of about 0.15 is “optimum” for the 4-stage 
scheme in the sense that it allows for the largest u for stability. A plot of the 
Fourier footprint at its maximum u = 1.75, along with the corresponding time 
stepping stability boundary using 77 = 0.15 is shown in figure 6.1. This exercise 
was also performed for 2, 3, and 5-stage schemes. Although not shown, for these 
cases the maximum CFL numbers for stability turn out to be v — 0.7, 1.2, and 
1.9, respectively. Hence increasing the number of stages yields increasing benefit 
in stability, but the difference in maximum allowable CFL number between 4 and 
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Figure 6.1: Stability Boundary of 4-Stage Time-Marching Scheme, and 

“Worst-Case” Fourier Footprint of First-Order Spatial-Dif- 
ferencing Scheme for v = 1.75 



-6 -4 -2 0 2 4 

Real 


Figure 6.2: Stability Boundary of 4-Stage Time-Marching Scheme, and 

“Worst-Case” Fourier Footprint of Second-Order Spatial- 
Differencing Scheme for v = 0.87 
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5 stages is small relative to the extra work involved in computing the extra stage. 
Therefore, in the present paper the 4-stage scheme is used for all explicit Euler 
computations. Note that since the Fourier footprint is highly dependent on Mach 
number, the maximum allowable v is actually higher than 1.75 for M lower than 
100. For example, at M = 3 the maximum CFL number is about 2.2, while at M 
= 1 it is about 2.5, according to this linearized stability analysis. 

For second-order spatial differencing (L defined by (6.8)), the Fourier trans- 
form of the right-hand side of (6.6) is: 

D(»< f) )| (-2cos< (i) + icos2C W + j) + 

i(|i(«W)) ( 2s mC( i > - i S in2C< i >) + |d(»«>)| (-2cos<«)+ ( 619 > 

icos2C«> + |) + <^(*.")) (““C 0 '’ - 5sin2(''>) . 

A “worst-case” Fourier footprint is again accomplished using the same vari- 
ables as for first order, only this time the maximum v turns out to be 0.87 for 
the 4-stage scheme with t) = 0.15. A plot of the Fourier footprint along with the 
stability boundary is given in figure 6.2. Again, at lower M the maximum CFL 
number is less restrictive. For M = 3 the maximum v is about 1.3, while at M = 
1 it is about 1.5. 

6.2 Implicit Time-Marching 

Ignoring the viscous terms, the implicit factored equation (3.15) can be writ- 
ten as 

/ 8F*\ 

= — +*„G*). 


I+JAtS, 


m 


A t U (n) 


( 6 . 20 ) 


S(-L) = - 
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Recall that the left-hand side Jacobian derivatives are written with first-order ac- 
curacy, resulting in block tridiagonals for each coordinate-direction sweep. Split- 
ting the Jacobian derivatives into + and — parts and expanding the right-hand 


side, ( 6 . 20 ) becomes 

[I + JAt(6-A+ + 6+Af)] [I + + «+ Bf )] A,U<”> 

= -^<(F; * +i , y - F-_ ;j + G* J+i - GJj.j), 


( 6 . 21 ) 


where the first-order operators are defined as: £*"(•) = (*)k+i — ( 0 * ^ k (') = 

(.) fc _ and k represents the index i or j. Approximate inviscid Jacobians 

are currently employed on the left-hand side for Aj*" and Bj^ , as defined in (3.20) 
and ( 3 . 22 ). 

The right-hand side of (6.21) is identical to the right-hand side of the explicit 
form (6.2), since $A 3 on £-faces equals F*, on 77 -faces equals G*, and 

J = 1/A. Hence, after linearizing the right-hand side as described in the last 
section, ( 6 . 21 ) can be written: 


MA,U (n) = -LU, 


( 6 . 22 ) 


where L is given by (6.7) or ( 6 . 8 ) and 

M = [I + JA<(*r Ai + + SfAr)} [I + + 6} Bf)] . (6.23) 


The Fourier symbols for the first-order difference operators in (6.23) are 



(6.24) 


Linearizing the A; and B; terms on the left-hand side of (6.22) about the cell 
center using Roe-averaged values, then taking the Fourier transform of the whole 
equation, the following expression for the amplification matrix g results: 


(S(M)}(g-I) = 3(-L), 


(6.25) 
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where 3(-L) is given by (6.16) for first-order right-hand side spatial accuracy or 
(6.19) for second-order, and 9f(M) is the Fourier transform of M, given by 


3(M) = < I + 






- cosC (,) )+ 


*ll + 




— (e[ j) ) (i- 
air 9 ’ v 


(6.26) 




COS 


Equation (6.25) can be rearranged to give the generalized eigenvalue problem: 


[3( M) + 3(-L)] x = g [9(M)] x, 


(6.27) 


where g is any of the complex eigenvalues, and [3(M) + ^(-L)] and [3( M)] are 
complex 4x4 matrices. 

The stability characteristics of the implicit equation are determined by cy- 
cling through 17 of each of the frequencies C (i) and C (i) from 0 to 2tt inclusive. 
(The combinations when both C (i) and C (j) equal either 0 or 2 tt are excluded since 
they yield eigenvalues of 1.0 automatically for a consistent scheme.) The gen- 
eralized eigenvalue problem is solved using a subroutine from the International 
Mathematics and Statistics Library (IMSL). The maximum eigenvalue, the aver- 
age eigenvalue, and the smoothing factor are determined as a function of the CFL 
number v. The maximum eigenvalue serves as an indication of stability, its value 
must remain at or below 1.0 for stability. The smoothing factor is the maximum 
eigenvalue when max(C (i \ < (i) ) lies between tt/ 2 and 3ir/2. It corresponds to the 
damping of high frequencies, and serves as an indication of the possible effective- 
ness of a multigrid method if it was applied to the scheme. The average eigenvalue 
simply gives a measure of the mean of all the eigenvalues over the whole frequency 


range. 
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Again, the stability is a function of a large number of independent variables, 
including the parameters z/, M, a, 0, 6' d , and 4j*\ so a thorough numerical analysis 
is difficult. For first-order spatial accuracy, the “worst case” parameters of all the 
variations tested were determined to be: M = 100, a = 90°, 6 d = 22.5 , 0 g — 0 , 
and 0 = 0 (although 0 = 1 gives a similar result here). As can be seen in 
figure 6.3 in a plot of the eigenvalue parameters as a function of CFL number, 
the maximum */ insuring stability is 2.5. For the different conditions of M=100, 
a = 45°, 9' d = 45°, 9 g = 0°, and 0 = 0, the scheme is stable only up to v = 0.05, 
the lowest limit found, as shown in figure 6.4(a). However, at these conditions the 
stability is a strong function of 0. When 0 = 0.05, the maximum * for stability 
increases to 2.5, as shown in figure 6.4(b). Higher 0 give even larger maximum 
allowable CFL numbers for these particular parameter's. As was mentioned in 
Chapter 5, 0 is limited to be greater than 0.05 as an aid to convergence, as 
determined empirically for solutions using explicit time- marching. Now it is clear 
that this limiting is necessary from the standpoint of stability for implicit-time 
marching as well. 

All other variations in parameters that were tested give stability limits that 
are either equally or less restrictive than those of figures 6.3 and 6.4. Hence, with 
0 limited to be greater than 0.05, the CFL limit for the first-order implicit scheme 
is about 2.5. As was true for explicit-time marching, the stability limit derived 
here for M=100 is somewhat overrestrictive for flows at lower Mach numbers. In 
practice, at reasonably low Mach numbers (less than about 3 or so) the stability 
limit appears to be about v — 4. 

When second-order spatial differencing on the right-hand side is employed, 
the CFL limit for stability is more restrictive. Figure 6.5 shows the stability plot 
for M = 100, a = 90°, 6' d = 22.5°, B g = 0°, and 0 = 0 (although 0 = 1 gives 
a similar result here). The maximum v = 1.4 for stability for this case. For the 
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Figure 6.5: Eigenvalues as Function of CFL Number for Implicit Scheme, 

Second-Order, M = 100, a = 90°, 0 ^ = 22.5 , = 0 , (3 = 0 
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conditions M=100, a = 45°, 6' d = 45°, 9 g = 0°, and ft = 0, the scheme is stable 
only up to the lowest limit of v = 0.01, as shown in figure 6.6(a). When (5 is 
limited to be greater than 0.05, then the stability limit for this case increases to 
about 0.3, as shown in figure 6.6(b). Hence, since other combinations of variables 
tested give similar or better stability limits, the maximum v for stability for the 
second-order implicit scheme is about 0.3. Again, at lower Mach numbers this 
stability restriction can be relaxed somewhat. In practice, at Mach numbers less 
than about 3 the stability limit is about u = 2. 

Naturally the question arises as to whether some left-hand side approximate 
Jacobians other than (3.20) and (3.22) can be devised to give better stability 
properties for the implicitly-advanced 5-wave model. It turns out that if the 
approximate Jacobians are taken as: 




( B 0; +i =5{(w) >+ r |6(# '' ,)l ^}’ 

where |D| is defined in (6.11), then the Fourier transform of the linearized left- 
hand side operator is given by 


I j+i 


(6.28) 


(6.29) 


3(M) = <1 + 




M^+w^)} 


D(^’ ) )|(l-cosC (i) ) + 


*<l + 


M4 0 ) + «(*$”)} 


Uh 


|D(^)|(1- (6.30) 


cosC W) ) + *(|^(^ i) ))«nC (,) J. 
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Figure 6.6: Eigenvalues as Function of CFL Number for Implicit Scheme, 

Second-Order, M = 100, o: = 45°, 6'^ = 45°, 6 g = 0 
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Numerical analysis of the resulting generalized eigenvalue problem shows this 
method to be unconditionally stable for first-order spatially accurate computa- 
tions. This linearized analytical result is confirmed in practice as well. First-order 
results have been seen to be stable at CFL numbers as high as 1000, although the 
optimum v for convergence generally occurs between about 5 and 10. 

Unfortunately, use of (6.28) and (6.29) in conjunction with a second-order 
right-hand side for the 5- wave model proves in practice to be even less stable than 
when the grid-aligned Jacobians (3.20) and (3.22) are used. An example of a 
stability analysis is shown in figure 6.7 for the same conditions as figure 6.5, only 
now the Fourier transform of the left-hand side (6.30) is used in place of (6.26). 
The maximum v for stability is only 0.6, compared with 1.4 for the grid-aligned 
approximate left-hand side. 

In order to attempt to improve the stability properties of the 5-wave model 
for second-order accurate computations using the approximately-factored implicit 
scheme, the method is reformulated using second-order left-hand side spatial ac- 
curacy. This makes it necessary to solve a block-pentadiagonal system rather than 
tri diagonals during each sweep. The procedure outlined by Barth [29] is followed 
to determine an appropriate form for the approximate left-hand side Jacobians 
for second-order accuracy. The end result is given here for the ith row of the 
block-banded matrix for the first sweep in solving the Euler equations: 


. . . , 0, (N x )i,(N 2 )i, ((N 3 )i + I/C JAt)) , (N 4 )i, (N 5 )i,0, . 


(6.31) 
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(N ‘ )<= [(w) i+ r5l 6(9 » ,) l^-il £,(# » i,) 

A similar form results for the second sweep. A stability analysis is carried out 


( 6 . 35 ) 

( 6 . 36 ) 


using 


M = [I + JAt(drA? + df A." )] [I + JAi^ + d/B; )] . 


( 6 . 37 ) 


The df are defined to be the following second-order difference operators: 


a^(-) = -|(-)* + 2 (‘)*+i _ 2^ k+2 

(•) = |(-)fc - 2 (-)*-i + 2(')*- 2 ’ 


( 6 . 38 ) 


where Jfc represents the index i or j. The Fourier symbols of these operators are: 

3 

+ 2e’' - ~e-' 

( 6 . 39 ) 

2 ' 2 

The Fourier transform of M, with terms linearized about the cell-center using 
Roe-averaged values, becomes: 


^ + )=-|+2e‘<-§e' 2t '' 


3(M) = jl + 
3 


{u;(dg^) + t*>(0fl^)} 


|D(^ ) )|(- 2c ° s C (i) + ^cos2C (i) 


+ 


*<l+ 


{ w (^) + 


Uh 


D(^ j) )| (-2cos£ (j) + -cos2C (j) + 2 ) + 


( 6 . 40 ) 


( 2sinC(i) ■ \ s[n2 ^ j) ) }• 


This linearized analysis shows the block pentadiagonal approximate-factorization 
method to be unconditionally stable in conjunction with a second-order spatially 
accurate right-hand side. In practice, however, the method is found to be stable 
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only up to CFL numbers of about 100 or so. This disagreement between theory 
and practice arises as a result of the highly nonlinear nature of the 5- wave model. 
The stability analysis is only a linearized approximation. Although the practical 
maximum allowable CFL number of 100 represents a dramatic improvement over 
the value of 2 that results when the grid-aligned left-hand side is employed, the 
optimum rate of convergence for steady-state calculations is still obtained for 
relatively low CFL numbers between about v = 2 and v = 6. 

Example convergence plots are shown in figure 6.8 for second-order spatially 
accurate computations of an inviscid shock reflection off a flat plate (M oo — 2.9) 
using the 5- wave model. When the grid-independent approximate Jacobians (pen- 
tadi agonal LHS) are employed, CFL numbers up to about 100 can be used, but 
the optimum rate of convergence is obtained at a CFL number of v — 3 for this 
case. Results with v = 2 are also shown for direct comparison with the result 
using the grid-aligned left-hand side Jacobians (tridiagonal LHS). Due to the in- 
consistency inherent in the use of grid-aligned left-hand side linearizations with 
the grid-independent right-hand side, the rate of convergence of the latter method 
tapers off as the solution approaches steady-state. However, both methods reduce 
the L 2 -norm of the residual to 10 -8 in about the same number of iterations. 

In summary, the stability properties of the 5-wave model, advanced implicitly 
in time using an approximate-factorization algorithm, can be improved dramati- 
cally through the proper choice of approximate Jacobians. However, for simplicity 
and consistency all implicit computations in the present study employ the grid- 
aligned approximate Jacobians described in Chapter 3. Although the correspond- 
ing maximum CFL numbers are more restrictive, they are usually high enough 
for most problems to allow for an acceptable rate of convergence. On the other 
hand, if the 5-wave model was to be used for higher Mach number flows where the 
stability restrictions using the grid-aligned left-hand side are more severe, or in 



log(L 2 -norm Res) 



Figure 6.7: Same as Figure 6.5, Except Approximate Jacobians (6.28) 

and (6.29) Used on Left-Hand Side 



Figure 6.8: Convergence Histories for Implicit A-F Scheme Applied to 

Second- Order Accurate Euler Computation Using 5- Wave 
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time-accurate computations where each mesh cell is advanced at a constant time 
step, it would probably be necessary to employ the grid-independent left-hand side 
terms in order to be able to advance the implicit solution at a reasonable rate. 
For second-order computations, this would entail the use of a block pentadiagonal 


solver. 



CHAPTER 7 


MONOTONICITY ANALYSIS 


The method for analyzing the monotonicity of the two-dimensional Euler 
equations is derived from considerations of the scalar convection equation Ut 4- 
au x = 0. The results of the Euler equations analysis axe considered to be valid for 
the Navier-Stokes equations as well, since the viscous terms add dissipation which 
tends to mitigate numerical oscillations that may occur near regions of strongly 
varying gradients. 

7.1 Scalar Wave Equation 


The one-dimensional scalar convection equation is written in finite-volume 
form, with forward- Euler time stepping ( i is a given cell bordered by (t — 1) to 
the left and (i + 1) to the right): 

(7.1) 

Here /i+ 1/2 and fi- 1/2 are the fluxes on the (i + 1/2) and (i — 1/2) cell faces, 
respectively, and Ax is the distance between the gridpoints. Consider now a com- 
putational stencil in which only depends on u ^ \ *4+i> and Godunov 

[3] showed that one way to insure that spurious oscillations do not develop is to 
require that variations in u in each cell in the computational stencil causes a vari- 
ation in the same direction in cell i. In other words, if increases, then u» 
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should also increase, or at worst remain unchanged. A similar requirement holds 
for changes in the (i + l)st cell. These requirements can be written 


du 


(n+l) 


Q ( n ) _ 

du t+1 


> 0 


and 


du\ 


(n+l) 


8ut\ 


> 0 


(7.2) 


or, since /t+ 1/2 is identical to f(ui+ i,Uj) and is identical to 

dfi-i 


dfax 
+ J < 0 


and 


> 0. 


(7.3) 


dui+i ~ du {- 1 

A third restriction is /du^ > 0, but this merely limits the time step. 

As an example of the usefulness of this analysis, consider first-order upwind 
differencing, which is already known to be monotone: 


fi+i = ~a( u »+i + «i) - ^M(*i+i - «i) 

/i_i = 2 °( u ‘ 1) — 2 l a l^ u * Ui ~ 1)' 


(7.4) 


Here, df i+1 / 2 /dui + 1 = |(a - |a|), which is non-positive, and = 

i(a + |ct|), which is non-negative. Hence first-order upwind differencing satisfies 
(7.3), as expected. A counter-example is central-differencing, which is already 
known not to be monotone: 


fi+l = 2 fl ( U *+l +“<) 

, 1 / x 

fi-\ = ~ a (Ut + 


(7.5) 


Here, it is seen that (7.3) can never be satisfied except for the degenerate case of 
a = 0, since both 5/i+i/2/^' u »+i an d ^ u i-\ — 2 a * 


7.2 Euler Equations 

The ideas from the last section extend in a straightforward manner to the 
one-dimensional Euler equations. These can be diagonalized yielding three non- 
linearly coupled convection equations, each which describes the convection of a 
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“characteristic variable.” Satisfying requirement (7.2) for each of these variables 
then implies monotonicity for the numerical scheme. 

Extension of this analysis technique for use with the two-dimensional Euler 
equations is not as simple because the equations are not diagonalizable in general, 
as discussed in Chapter 2. However, the influence of variations in the conserved 
variables can be decoupled locally, as described below. First, the two-dimensional 
Euler equations are written in finite-volume form with forward-Euler time step- 
ping, and it is assumed for simplicity that the mesh is made up of square cells 
with face length As: 


jj(n+l) Tj( n ) 


At 

As 




*+ j ,j 




+ * 


*.i+ 2 



(7.6) 


U are the conserved variables and $ are the fluxes per unit face-length on the cell 
faces, given by (5.10). It is assumed that the computational stencil is made up of 
only (», j), (i + 1 , j), (» - l,j), (i,j + 1), and (i, j - 1), so that in one time step 
U t - • is only a function of its initial value and the values in the four immediate 
neighboring cells. Thus this analysis only applies to a spatially first-order accurate 
scheme. 

Now, instead of one equation, there are four coupled equations, and the quan- 
tity for example, is not a single variable but a 4 x 4 matrix. 

The matrices 0u£ +1) /dU ( * n) , where k = (z + 1, j), (* - 1, j), (*, j + 1), (*,J ~ 1), 
are termed the “influence matrices.” The four eigenvalues of each influence ma- 
trix represent the change of certain characteristic variables at (i ? j) caused by 
the change in those same four variables in the corresponding neighboring cell. 
Since the equations describe the convection of these characteristic variables, non- 
negativity of the eigenvalues of each influence matrix implies monotonicity for the 
corresponding set of characteristic variables. Note, however, that the characteris- 
tic variables associated with the cells (z ± 1 , J? ) are not the same as those for cells 
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(i, j ± 1), except for the entropy; therefore it is not clear what the combined effect 
of these monotonicity properties is. Nonetheless, this approach is utilized to help 
define a limiting procedure for reducing oscillations in two-dimensional solutions. 
The conditions equivalent to (7.2) for the two-dimensional Euler equations 


are 


e.v. | ^5-c-^ > 0 k = (i + - 1. X ) ( 7 - 7 ) 

aul n) } 




or, equivalently, 


e.v. 


a# 




dV i+1<j 


< o 


te) - 0 


(7.8) 


where e.v.(-) represents “the eigenvalues of (•)”• If the gnd-normal angle 9 g is 
varied over the full range of possible angles, then satisfying all four inequalities 
in (7.8) is redundant. Satisfying the two inequalities on opposing faces (say the 
(i+1 12, j) face and the (z— 1/2, j) face) is then sufficient to insure this monotonicity 


property. 

In order to proceed with the monotonicity analysis for the 5-wave 
independent model, (5.10) is written in slightly different form: 






\ + *«) ■ -\ { [R1 1 (A-]c OS («i - «,)|n*+ 


grid- 


(7.9) 


[R*] is the matrix of wave vectors acting in the ^-direction (Ri. 2,4,5 in (5.11)). 
[R**] is the corresponding matrix of wave vectors acting in the (9' d + Indirection 
(only the shear vector, R 3 in (5.11), is used in the 5-wave model). Then [A*] = 
diag {q' d + a, q' d - a,q' d ,q' d ) and [A**] = diag(f^ + a, f' d - a,r' d ,r' d ) are the corre- 


sponding wavespeeds in those same directions. The wavestrengths are 

r££+ fi-fr ( Aucosflj + Arising) ' 

6* - A “ (Ancos 6' d + Ausin Q' d ) 

° - ±(a*Ap-Ap) 

| (-Ausin 9 d -I- Avcos 9 d ) 


(7.10) 
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n** = 


o 

o 

0 


(7.11) 


LOS - 1)| (Atxcos^ + Avsin0^) . 

The present study concentrates on the variations in O* and O* . It is assumed 
that the wave vectors, wavespeeds, and 8' d are constant, and all variables are taken 
as the Roe-averaged values. In addition, (3, although a function of A u, Au, and 
A p, is considered constant in order to simplify the final expressions. With these 
assumptions, this becomes a linearized analysis and one can obtain 

50* 




au R 


fin** 1 


(7.12) 


In a similar fashion, for the (i — 1/2) face 


an* 


au L 


an** l 


(7.13) 


The monotonicity constraints are 



|<o 

au i+1 > 


| > o. 

9Ui_! > 


(7.14) 


The derivative matrices for the 5- wave model can be found as dCl /$Ur 

( 7 - l)u , 0CO3 e’ d (l-l)v , 7^1 

oaV I oA 2a 2 ' 2d 2d 3 


r S^(i 2 +* ! )-# 

^(i 2 +« 2 )+# 

_ ( T -l)(u 3 +t> 3 ) 

1 2d 3 

r'j 


an* 


2d 3 1 2d 

(7-l)t* _ /3cosO^ 
2d 3 2d 

d 3 

sin St 


2d 3 1 2d 

2d 3 2d 

d 3 

cos St 


2d- 

2=1 

2d 3 

_il =12 

a 2 

o 


, (7.15) 


dU R 


0 0 0 0 

0 0 0 0 

0 0 0 0 

(j9-l)cosg^ (/j-i^sinS^ 0 J 


(7.16) 
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and (da*/du L ) = -(an*/au a ), and (an**/au L ) = -(dn**/d u R ). 

The monotonicity analysis is carried out numerically. The Mach number M, 
flow angle a, and 0 are chosen, then 8 g and 0' d are each varied independently be- 
tween -90° and 90° with incremental changes of tt/ 32. Eigenvalues are computed 
for each condition. They are usually real numbers, but can also be complex con- 
jugates; in such cases only the real parts are considered. If they meet the criteria 
of equation (7.14), then monotonicity is preserved at that condition. It turns out 
that plotting (0’ d - a) vs. (8 g - a) removes the dependence on a (in other words, 
plots are the same regardless of the value of a). 

A sample plot is shown in figure 7.1(a). The conditions are M=3, ^=0.95. 

There are two very small regions where monotonicity is preserved. (Note that 

some points may be missing from these monotonicity plots wherever the eigenvalue 

solver does not converge within a specified number of iterations. However, we are 

more interested in general regions than in specific points.) As a specific example, 

from the figure it is seen that the scheme is monotone for approximately 30 < 

0'-a< 75° when 9 a - a = 75°. These example allowable conditions are sketched 
d y 

in figure 7.1(b). 

It is also evident from figure 7.1(a) that if ( & g at) lies between roughly 60 
and 60°, then no 8' d chosen will insure monotonicity. Other 0's less than 1.0 
produce similar plots. Only when 0 = 1.0 is there always some 6' d that will yield 
a monotone scheme, as shown in figures 7.2(a) through (c). Here, the diagonal 
where 8 d = 6 g corresponds to the grid-aligned scheme. The effect of Mach number 
is also shown in the figures. At low Mach numbers, only the grid-aligned method 
is monotone (t.e. 6' d must = 8 g and 0 must = 1.0), while at higher supersonic 
Mach numbers the monotone region is extended slightly from the M=3 case. 

It is clear from this analysis that the restrictions on allowable 0' d for a mono- 
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100 


e g -a, deg 


a) Monotonicity Regions 



b) Approximate Allowable d' d When 6 g — a = 75° 
Figure 7.1: Monotonicity for M = 3, f3 = 0.95 
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tone scheme given by this analysis are quite severe, if not impossible to meet. For- 
tunately, in practice it appears that the restrictions on 9' d can be relaxed somewhat 
while still maintaining reasonably non-oscillatory solutions near discontinuities for 
a wide variety of flows. 

Through an extensive amount of numerical experimentation with actual solu- 
tions to the Euler equations, the following observations have been made regarding 
reducing the oscillatory behavior of the grid-independent model to an “acceptable” 
level: 

(1) When M > 1, best results are obtained when 9' d is limited to lie between 
a + K * sign(0 s — a) and 0 g . K is a small number for lower Mach numbers and is 
larger for higher Mach numbers. 

(2) When M « 1, 0' d does not need to be restricted, except in a very small 
region (see (3) below). Between M = 0 and M = 1, the allowable region is 
transitioned smoothly between the subsonic and supersonic cases. 

(3) In the boundary layer region of Navier-Stokes solutions, odd-even point 
decoupling can occur when 6' d is taken as (0 a i -j), and a » 0' d . This condition 
occurs on grid interfaces in the boundary layer that are aligned with the flow 
direction, and is due to the fact that all components of the 0^-wavespeeds in 
the grid direction equal zero, and the ( 6 d + shear wave has an extremely small 
wavespeed. Hence the dissipation is very small, and the result is essentially central- 
differencing in that direction. By limiting the angle 6' d to lie outside of a small 
region near (0^ — c*) = 0 at (0 9 — a) = i90°, this decoupling can be alleviated. 
Numerical examples of viscous flows both with and without 0^-limiting will be 
given in Chapter 8. 

An attempt has been made to parameterize the “monotonicity regions’ in 
accordance with the three observations made above. The empirically-generated 
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0 j-limited regions for four different Mach numbers are shown in figures 7.3(a) 
through (d). It should be stressed that the determination of these regions is 
based only loosely on theory and primarily on numerical experimentation. The 
following empirical scheme has been found to give good results for a wide variety 
of problems. It is by no means deemed to be the best scheme for improving the 
monotonicity properties of the 5-wave model. First, some variables are defined: 

(7.17) 

(7.18) 

(7.19) 

(7.20) 

(7.21) 

(7.22) 
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y 7 : -K /max 


(f) ! -{<*»-“> + if ’° 


i/s = -y? 


V9 = + 


yLax [(£)’- {(« s - a) -|} ,1 


2/io = -2/9 

2/h = min(t/ 6 ,max(2/4,y9)) 

2/12 = max ( 2 / 5 , min( 2 / 3 , ys )) • 

The allowable regions are then taken as: 

(9' d - a) > max(yi,i/ 7 ) 

(9' d - a) < min(j/ 2 >yio) 

yu < (#d - «) < yi2- 


(7.23) 
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Figure 7.3: Empirically-Determined Mon< 
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If ( 0' d - a) does not lie within one of the allowable regions, then it is limited 
to either y u or y 12 , whichever is closer. 



CHAPTER 8 


TWO-DIMENSIONAL RESULTS 


8.1 Euler Computations 

Both first-order and second-order two-dimensional Euler computations are 
performed for several different cases in this section. Unless otherwise noted, 
second-order computations do not utilize any type of limiting of higher-order terms 
for either the grid-aligned or grid-independent computations. This is done in order 
to avoid confusion with and separate the effects of the ^-limiting procedure. 

8.1.1 Shock Reflection 

The M = 2.9 inviscid shock reflection case is computed on a Cartesian mesh 
4.8 units wide by 1.6 units high. An oblique shock enters the domain from the 
upper left corner, reflects off the bottom wall, and exits out the right end. The 
flow is turned through an angle of 10° by the incident shock. The nondimensional 
boundary conditions (nondimensionalized by combinations of p^ and doo) are: 
at inflow p = 1, pu = 2.9, pv = 0, and pE = 5.9907; at the top boundary 
p = 1.6328, pu = 4.3272, pv = -0.7630, and pE = 9.5091; at the back boundary, 
outflow conditions are set by second-order extrapolation from the interior; at 
the body along the bottom wall, simple reflection boundary conditions are used. 
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Computations are performed using explicit four-stage time-marching. They are 
initiated from freestream conditions and are run until the L 2 -norm of the residual 
of all four equations drops below 1 x 1(T 12 . Computations are performed on two 
different mesh sizes: a 49 x 17 mesh and a 97 x 33 mesh. 

First-order computations using the grid-aligned solver, run at a CFL number 
of i/ = 2.2, converge in 158 iterations on the coarse mesh and 258 iterations on the 
fine mesh. Nondimensional pressure contours and pressure values along three j = 
constant cuts (left to right through the mesh) are shown in figures 8.1 and 8.2. 
First-order results using the 5-wave model are also obtained using i/ = 2.2. The 6 d 
values are frozen using the following procedure: they are computed every iteration 
for the first 20 iterations, then only once every 20 iterations until the log of the L 2 - 
norm of the residual drops to below -3.5. After this, the 6 d values remain frozen. 
The 5-wave model is run on the coarser mesh both with and without the e' d - 
limiting derived in Chapter 7. The solution converges in in 245 iterations {6' d not 
limited) and 191 iterations (0' -limited). On the finer mesh the solution converges 
in 319 iterations (^-limited). The results are given in figures 8.3 through 8.5. On 
both meshes, the 5-wave model gives sharper shock wave resolution than the gnd- 
aligned scheme. The ^-unlimited method produces the sharpest resolution, but 
at the cost of oscillatory behavior near the discontinuities. When the ^-limiting 
procedure is employed, 5-wave model results are still significantly sharper than 
the grid-aligned results, and appear to be monotone as well for this problem. 

The grid-aligned model on the 97 x 33 grid converges in about 14.1 CPU 
seconds on the Cray 2 computer, while the 5-wave model requires about 25.4 
CPU seconds (^-limited). Thus, the 5-wave model is about 1.8 times more costly 
than the grid-aligned model to reach the same level of convergence for this case. 
Per iteration, the 5-wave model is about 1.5 times more costly. A significant 
percentage of this cost is associated with the limiting of the 0' -directions. Without 
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b) Pressure Along j = Constant Lines 

Figure 8.3: Shock Reflection, First-Order, 5- Wave (^-Unlimited), 

49 x 17 



pr essure 


a) Pressure Contours 



b) Pressure Along j = Constant Lines 


Figure 8.4: Shock Reflection, First-Order, 5- Wave (^-Limited), 49 x 17 
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^-limiting, the 5-wave model is only about 1.2 times more costly than the grid- 
aligned model per iteration. 

Second-order results for this problem are obtained using a CFL number of 
v = 1.1. Grid-aligned model results are shown in figures 8.6 and 8.7 on the two 
grids. Since no limiting of higher-order terms is employed, there are undershoots 
present in the solutions forward of each shock wave. These solutions converge 
the L/ 2 -norm of the residual of all four equations to 1 x 10 12 in 249 and 423 
iterations, respectively. Second-order 5-wave model results are given in figures 
8.8 through 8.10. A 0^-unlimited solution on the 49 x 17 mesh (figure 8.8) shows 
very sharp shocks, but a large amount of oscillatory behavior is also present aft 
of the shocks. This solution converges in 494 iterations. When ^-limiting is 
used, these oscillations are reduced in magnitude, while still maintaining slightly 
sharper shock resolution than the grid-aligned method as shown in figure 8.9. In 
the figure, the wiggles in the “3.6”-level contour are caused by oscillations which 
are still present downstream of the reflected shock. The 5-wave d^-limited solution 
on the 97 x 33 grid is shown in figure 8.10. Again, results are slightly sharper than 
the grid-aligned method on the same grid and additional oscillations are visible 
downstream of the reflected shock. The ^-limited solutions converge in 391 and 
659 iterations, respectively, on the coarse and fine grids. 

For the shock reflection problem, the 5-wave model appears to be a viable al- 
ternative to the grid-aligned model for first-order computations. Wfien ^-limiting 
is employed, results appear to be free from spurious oscillations, and shocks are 
captured with fewer interior points than the grid-aligned method. In fact, the 
49 x 17 5- wave solution gives comparable resolution to the grid-aligned solution 
on a mesh with four times as many mesh cells. For second-order computations, 
very little advantage of the 5-wave model over the grid-aligned model is seen. 
Results are slightly sharper, but the extra cost for the method may outweigh its 
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b) Pressure Along j — Constant Lines 
Figure 8.7: Shock Reflection, Second-Order, Grid-aligned, 97 x 33 








pressure 


a) Pressure Contours 



b) Pressure Along j = Constant Lines 


Figure 8.10: Shock Reflection, Second-Order, 5- Wave (^-Limited) 

97 x 33 
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small advantage. Also, small oscillations still plague the solutions. 


8.1.2 Ramp Flow in a Channel 


Euler computations are performed for supersonic channel flow with a 15° 
finite-length ramp at an inflow Mach number of 2.0 on a 49 x 17 grid. This case 
was first used to test grid-independent flow solvers by Levy et al. [24]. The grid 
is 3 units long by 1 unit high. The ramp begins at x = 0.5 units and is 0.5 units 
long. At the end of the ramp, the bottom channel wall again becomes parallel 
to the top wall. A picture of the grid is given in figure 8.11. The boundary 
conditions are uniform freestream inflow at the left, second-order extrapolation at 
the right, and simple reflection conditions at the top and bottom walls. For this 
case, some extra dissipation necessary to prevent oscillations in the 5-wave model 
solution near the bottom wall in the region of the expansion fan was provided by 
using the grid-aligned scheme in combination with the simple reflection boundary 
conditions to obtain the flux at the lower wall. 

A first-order computation using the grid-aligned method is shown in figure 
8.12. As in the shock reflection case, an explicit four-stage time-marching scheme 
is employed using a CFL number of i/ = 2.2. Figure 8.12(a) shows Mach number 
contours, while figure 8.12(b) shows Mach number values along two j = constant 
lines (left to right through the mesh). It can be seen that this grid-aligned first- 
order result smears the primary shock a good deal and barely shows evidence 
of the two reflected shocks. The 5-wave model with no 6' d limiting, shown in 
figure 8.13, shows extremely sharp resolution of the shocks. However, significant 
oscillations are also present in the solution. When ^-limiting is employed, the 
oscillations essentially disappear, but much of the resolution of the 5-wave model 
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Figure 8.11: Ramp Flow 49 x 17 Grid 




b) Mach Number Along j = Constant Lines 


Figure 8.12: Ramp Flow, First-Order, Grid-aligned 
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b) Mach Number Along j = Constant Lines 
Figure 8.13: Ramp Flow, First-Order, 5-Wave (^-Unlimited) 




X 

b) Mach Number Along j = Constant Lines 
Figure 8.14: Ramp Flow, First-Order, 5- Wave (^-Limited) 


110 


is also lost. These results are shown in figure 8.14. They are still sharper than 
the grid-aligned results. The results from figures 8.12, 8.13 and 8.14 converge the 
L 2 -norm of the residual of all four equations to 1 x 10“ 12 in 263, 346, and 315 
iterations, respectively. The 5-wave model uses the same strategy for freezing 6<i 
described for the shock reflection case. 

Second-order solutions are obtained with u = 1.1. Results using the grid- 
aligned model are shown in figure 8.15, while 5-wave model results without and 
with ^-limiting are shown in figures 8.16 and 8.17, respectively. When no Q d - 
limiting is employed, results are very sharp but oscillatory. The ^-limited solu- 
tion is significantly less oscillatory, and is slightly sharper than the grid-aligned 
solution. However, as for the second-order shock reflection case, the amount of 
benefit probably does not outweigh the disadvantages in this case. 

8.1.3 Oblique Supersonic Shear 


The oblique shear wave case is computed on a 61 x 21 Cartesian mesh 3 units 
wide by 1 unit high. Fluid enters the domain from the lower face at a 45° angle 
and exits out the top face. To the left of the shear wave M = 1.812, while to the 
right M = 1.510; there is one transition cell where M = 1.661. There is no pressure 
change through the shear. The nondimensional boundary conditions are: to the 
left of the shear p — 1, pu — 1.282, pv = 1.282, and pE = 3.427; to the right of 
the shear p — 1, pu — 1.068, pv = 1.068, and pE — 2.926; at the transition zone 
p = 1, pu = 1.175, pv = 1.175, and pE = 3.165; on the top and right boundaries, 
outflow conditions are obtained using second-order extrapolation. 

The exact solution is shown in figure 8.18. Figure 8.18(a) shows Mach number 
contours while figure 8.18(b) shows Mach number values along three j = constant 
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b) Mach Number Along j — Constant Lines 
Figure 8.15: Ramp Flow, Second-Order, Grid-aligned 



a) Mach Number Contours 





b) Mach Number Along j — Constant Lines 
Figure 8.16: Ramp Flow, Second-Order, 5- Wave (^-Unlimited) 
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b) Mach Number Along j = Constant Lines 
Figure 8.17: Ramp Flow, Secoud-Order, 5- Wave (^-Limited) 



b) Mach Number Along j = Constant Lines 
Figure 8.18: Oblique Shear, Exact Solution 
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lines through the mesh. All computations are performed using four-stage explicit 
time-marching. When initiated from the exact solution, a first-order Euler com- 
putation using the grid-aligned method smears the shear wave significantly, as 
shown in figure 8.19. The first-order 5-wave model, with ^-limiting, converges 
to the results shown in figure 8.20. Results are significantly sharper, although 
there is still some oscillatory behavior near the discontinuity. For this particular 
example, the initial condition for the 5-wave model is the grid-aligned solution. 
The Q& values are recomputed once every 30 iterations until the log of the L 2 -norm 
of the residual drops below -3.5, after which they remain frozen. It should be 
noted that if /3 is allowed to be zero (instead of being restricted to be greater than 
0.05) and the outflow boundary conditions are extrapolated from 45 upstream 
(instead of in a grid-aligned fashion), then the 5-wave model can maintain the 
exact solution in one iteration when the exact solution is the initial condition. 

Second-order solutions using the grid-aligned and 5- wave models are shown in 
figures 8.21 and 8.22. The grid-aligned results are now fairly sharp, comparable in 
width to the first-order 5-wave results. The second-order 5-wave results are even 
sharper, however. The shear wave is now resolved with almost no spreading at 
all. Because this is a fairly weak shear wave case, the grid-aligned model does not 
yield any noteworthy pressure distortions as a result of misinterpreting the oblique 
shear wave. Both the grid-aligned and 5-wave models produce pressure fields in 
error from the exact solution by less than 0.5%. This 45° supersonic shear wave 
case demonstrates an advantage of the 5-wave model over the grid-aligned model. 
For both first and second-order computations, the model resolves the oblique shear 
wave with significantly fewer interior points than the grid-aligned method. 
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b) Mach Number Along j = Constant Lines 
Figure 8.19: Oblique Shear, First-Order, Grid-Aligned 
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b) Mach Number Along j = Constant Lines 


Figure 8.20: Oblique Shear, First-Order, 5- Wave 
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b) Mach Number Along j = Constant Lines 
Figure 8.21: Oblique Shear, Second-Order, Grid-Aligned 
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b) Mach Number Along j = Constant Lines 
Figure 8.22: Oblique Shear, Second-Order, 5- Wave 
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8,1.4 Supersonic Flow Over an Airfoil 


Euler computations are performed for the NACA 0012 airfoil at M 1.2, 
a = 0° on several different grids. The finest is a 257 x 73 O-mesh with an outer 
boundary extent of 20 chords and an average minimum spacing at the body of 
0.0031 chords. It is shown in figure 8.23. Two coarser meshes, a 129 x 37 and 
a 65 x 19, result from removing every other point from the next finest mesh. 
The inflow /out flow boundary conditions in the farfield for supersonic flow are 
simply freestream values when the grid boundary cell is an inflow boundary (t.e. 
the outward-pointing grid-face normal points into the freestream direction), and 
are extrapolated from the interior using second-order interpolation when it is an 
outflow boundary. The boundary conditions at the body are simple reflection 
conditions, and periodic boundary conditions are enforced where the grid meets 
itself behind the airfoil. 

These airfoil results, as well as all results to follow in the remainder of the 
paper, are computed using the implicit approximate-factorization time-stepping 
procedure. The CFL numbers at which the solutions are advanced are taken in 
accordance with the analysis performed in Chapter 6. Results are not converged 
to machine zero in general; sufficient convergence is assumed when the L 2 -norm 
of the residual drops by at least 4 orders of magnitude and/or the lift and drag 
values settle down and do not vary significantly with further iterations. 

At the conditions enumerated above, the NACA 0012 airfoil has a flowfield 
with a curved bow shock located in front of the airfoil and oblique shocks ema- 
nating from the trailing edge. Figures 8.24 through 8.26 show results for the first- 
order accurate grid-aligned model on the three successively finer grids. Shown are 
nondimensional pressure contours, plotted in increments of 0.05 (the freestream 



Figure 8.23: NACA 0012 257 x 73 O-Grid 




Figure 8.24: Supersonic Airfoil-Flow, First-Order, Grid- Aligned, 65 x 19 
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value of 1.0 is not plotted), as well as pressure along a j = constant grid line 
over the upper half of the airfoil. The j = constant grid line is taken at j = 11, 
j = 21, and j = 41 for the three grids. Its location is shown in figure 8.27. The 
two finer meshes produce reasonable-looking results, although the shock waves 
are very smeared. The coarsest mesh (65 x 19) produces very poor results, but 
this grid is so coarse that a good solution would not be expected in any case. 
First-order 5-wave model results (with 0^-limiting) are computed by restarting 
the corresponding grid-aligned solution with 9 d frozen. Results are shown in fig- 
ures 8.28 through 8.30. The shock waves are captured much more sharply by this 
method, although there are still some small oscillations present near the discon- 
tinuities. The 65 x 19 mesh is again too coarse to produce a reasonable-looking 

solution. 

In an effort to explore the grid-sensitivity of the grid-aligned and 5-wave 
models, computed drag values are plotted for each of the grids in figure 8.31. For 
a first-order scheme, it is expected that the results behave in a linear fashion when 
plotted vs. the inverse of the square root of the total number of gridpoints. This 
indeed appears to be the case for both methods (on grid sizes of 129 x 37 or finer). 
The drag coefficient is approaching about 0.0955 as the grid approaches infinite 
refinement. The 5-wave model produces drag values in closer agreement with this 
“correct” value on all three grid sizes tested. For example, on the 257 x 73 gnd, 
the 5-wave model gives a drag value 2.0% in error from the extrapolated “exact” 
value, while the grid-aligned result is 6.5% in error. 

Second-order computations using the grid-aligned model on the three gnds 
are given in figures 8.32 through 8.34. Results now have much sharper resolution. 
Since no limiting of higher-order terms is performed in these computations, some 
undershoots and overshoots are present near the computed shock waves. Second- 
order 5-wave model results are given in figures 8.35 through 8.37. These results are 



a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.29: Supersonic Airfoil-Flow, First-Order, 5-Wave, 129 x 37 
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a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.30: Supersonic Airfoil-Flow, First-Order, 5- Wave, 257 x 73 



Figure 8.31: Grid-Convergence Study for Supersonic Airfoil-Flow, First- 

Order 
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a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.32: Supersonic Airfoil- Flow, Second-Order, Grid- Aligned, 65 x 19 


a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.33: Supersonic Airfoil-Flow, Second-Order, Grid-Aligned, 

129 x 37 



Z) 

cn , , 
in 1 • 1 
<u 


.7 1 

-1.5 



x/c 


a) Pressure Contours b) Pressure Along j - Constant Line 

Figure 8.34: Supersonic Airfoil-Flow, Second-Order, Grid- Aligned, 

257 x 73 


a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.35: Supersonic Airfoil-Flow, Second-Order, 5- Wave, 65 x 19 
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a) Pressure Contours b) Pressure Along j = Constant Line 

Figure 8.36: Supersonic Airfoil-Flow, Second-Order, 5- Wave, 129 x 37 


a) Pressure Contours b) Pressure Along j - Constant Line 

Figure 8.37: Supersonic Airfoil-Flow, Second-Order, 5- Wave, 257 x 73 
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Figure 8.38: Grid-Convergence Study for Supersonic Airfoil-Flow. 

Second- Order 


126 


still somewhat sharper than the grid-aligned results, but there are also significant 
oscillations near the shocks. Computed drag coefficient values are plotted as a 
function of the inverse of the grid size in figure 8.38. Second-order results behave 
linearly when plotted this way (for a sufficiently fine mesh). The drag coefficient 
again extrapolates to about 0.0955 on an infinitely fine mesh for both methods. 
The 5-wave model gives better predictions for the drag on all three grids. On 
the finest, it deviates from the extrapolated value by 0.05% while the gnd-aligned 
method deviates by 0.22%. 

In summary, use of the 5- wave model for supersonic flow over an airfoil results 
in sharper shock resolution and better airfoil drag prediction than the grid-aligned 
model. The increase in shock resolution is more dramatic for first-order computa- 
tional accuracy, although even second-order results show a noticeable difference. 
In spite of the fact that ^-limiting was employed, both first and second-order 
5-wave model results are somewhat oscillatory near the computed shock waves. 

The effect of limiting the higher-order terms in the second-order solution using 
the 5-wave model is investigated at this point. A computation is performed on 
the finest mesh with the “min-mod” limiter employed. This limiter is described 
in more detail by Anderson et al. [30]. Results are shown in figure 8.39. The 
shock waves are resolved with the same number of interior points as figure 8.37, 
but now the amount of oscillations near the discontinuities is reduced. Hence 
the use of standard limiters for second-order computations seems to be viable 
in conjunction with the 5-wave ^-limited model. While the limiter does not 
guarantee oscillation-free solutions for this model, it does seem to reduce the 
amount of oscillations present. 
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8.1 ,5 Subsonic Flow Over an Airfoil 


A grid convergence study is performed using the Euler equations for subsonic 
flow over an airfoil in order to investigate the effect of the 5-wave model on an 
airfoil flowfield where no shock waves are present. The conditions are M = 0.3 
and a = 1°, and results are computed on the same three O-meshes used in the 
supersonic airfoil flow study. Boundary conditions on the body are again simple 
reflection conditions, and periodic boundary conditions are applied where the mesh 
meets itself behind the airfoil. The inflow/outflow boundary conditions applied 
at the outer boundary for subsonic flow are obtained using characteristic theory. 
To summarize, the local Riemann invariants R+ and R" are constructed as 


where q g indicates the velocity normal to the outer grid boundary face and the 
superscripts “(*)” and “(o)” indicate that conditions are taken from inside and 
outside (= freestream) the grid boundary, respectively. The normal velocity and 
speed of sound at the face are then taken as 


(9g)f = ^( R+ + R ) 

o, = ^-i(R + -R-)- 


( 8 . 2 ) 


If > 0, the entropy „ and the tangential velocity (r,), are extrapolated from 

inside the grid boundary. Otherwise, they are taken as the freestream values. The 
density boundary condition is then obtained using 


Pi - 



(8.3) 


The u and v velocity components are constructed from (q g ){ and (r g ) f, and the 
energy is calculated using the equation of state. 
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Figure 8.39: Supersonic Airfoil-Flow, Second-Order, 5- Wave, 257 x 73, 

Min- Mod Limiter Employed 



Figure 8.40: Grid-Convergence Study for Subsonic Airfoil-Flow, First- 

Order 
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Figure 8.40 is a plot of computed drag coefficient vs. the inverse of the square 
root of the grid density for the grid-aligned scheme and the 5- wave model on the 
three grids using first-order spatial discretization. The “exact” Euler solution 
should give zero drag. The 5-wave model, restarted from the grid-aligned solution 
with frozen, gives a far better prediction of the drag than the grid-aligned 
scheme for all three grids. Entropy contours (where entropy is defined as p/(p 7 ) - 
1) for both methods are given in figures 8.41 and 8.42. Contour values plotted are 
in increments of 0.001 for both figures. These figures indicate significantly lower 
entropy production over the airfoil surface using the 5- wave model. The maximum 
entropy values for the grid-aligned model are 0.0303, 0.0251, and 0.0183 for the 
coarsest through finest meshes, respectively. For the 5-wave model the maximum 
values are only 0.0052, 0.0028, and 0.0018. 

It is believed that the difference in entropy levels is due to the different ways 
that the two models interpret the flow near the stagnation point of the airfoil. 
Near the stagnation point, the flow undergoes very rapid turning with relatively 
small changes in pressure. The grid-aligned model interprets this turning to be 
in part due to the action of + and - acoustic waves with nearly offsetting Ap’s. 
Because the local flow is subsonic, the wavespeeds associated with each of these 
acoustic waves are of opposite sign, so the flux computed at interfaces near the 
leading edge is assigned a pressure which is too high or low by an amount A p. This 
results in increased entropy generation. In contrast, the 5-wave model interprets 
the rapid turning near the stagnation point to be due primarily to the action of a 
(0^ _|_ 2) shear wave, which has no associated pressure jump across it. Numerical 
entropy generation is lower as a consequence. 

Total pressure loss values (= (ptoo — Pt)/ptoo) along the upper surface of the 
airfoil are plotted for the grid-aligned and the 5-wave models in figures 8.43(a) 
and (b). The total pressure loss approaches zero for both models as the grid is 
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b) 129x37 



c) 257 x 73 



Figure 8.41: Subsonic Airfoil- Flow, Entropy Contours, First-Order, Grid 

Aligned 


a) 65 x 19 



b) 129 x 37 



c) 257 x 73 



Figure 8.42: Subsonic Airfoil-Flow, Entropy Contours, First-Order, 

5- Wave 
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Figure 8.43: Subsonic Airfoil-Flow, Upper Surface Total Pressure Loss 

First- Order 
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refined, as expected. Near the leading edge the grid-aligned model shows much 
higher magnitudes than the 5-wave model, reflecting the higher entropy generated 
there. Over the majority of the surface of the airfoil on the finest mesh, however, 
the grid-aligned model shows lower total pressure loss values. It appears that 
these are approaching zero at a faster rate than those from the 5-wave model as 
the grid is refined. The reason for this is unknown. 

Drag coefficient values from the grid-sensitivity study for second-order accu- 
rate computations are plotted in figure 8.44. Again, the 5-wave model gives lower 
values than the grid-aligned model, but in this case the difference is not so dra- 
matic as for first-order. This is reflected in the entropy contours as well, shown 
in figure 8.45 for grid-aligned and 8.46 for the 5-wave model. (Contour values 
plotted are in increments of 5 x 10" 5 for both figures.) The 5- wave model appears 
to have only slightly lower entropy values overall. The maximum values for the 
grid-aligned method are 0.0018, 0.0011, and 0.0006 on the coarsest to finest grids, 
respectively, while the maximum values for the 5-wave model are 0.0015, 0.0009, 
and 0.0007. Total pressure loss values on the airfoil upper surface are plotted in 
figures 8.47(a) and (b). In general, results are similar for the two methods. The 
5-wave model gives somewhat higher values just aft of the leading edge than the 
grid-aligned model on all three grids, but values near the trailing edge are lower. 
For completeness, the pressure contours resulting from the second-order fine mesh 
5-wave model solution are shown in figure 8.48. Pressure levels in increments of 
0.005 are plotted. The grid-aligned result on the finest mesh is identical to this 
figure, so it is not shown. 

In summary, first-order computations using the 5- wave model give lower nu- 
merical entropy generation in a subsonic Euler computation over an airfoil. As 
a result, the drag prediction is nearer to the correct result of zero. Second-order 
computations using the 5-wave model still give better drag prediction than the 
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Figure 8.44: Grid-Convergence Study for Subsonic Airfoil-Flow, Second- 

Order 
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b) 129 x 37 




c) 257 x 73 



Figure 8.45: Subsonic Airfoil-Flow, Entropy Contours, Second-Order, 

Grid- Aligned 
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a) 65 x 19 


b) 129 x 37 




c) 257 x 73 



Figure 8.46: Subsonic Airfoil-Flow, Entropy Contours, Second-Order, 

5- Wave 
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Figure 8.47: Subsonic Airfoil-Flow, Upper Surface Total Pressure Loss, 

Second-Order 
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Figure 8.47: Concluded 



Figure 8.48: Subsonic Airfoil-Flow, Pressure Contours, Second-Order, 

257 x 73 
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grid-aligned method, but the difference is not as great as for first-order. Entropy 
levels are similar for the two methods. Hence there does not appear to be a sig- 
nificant advantage to using the 5-wave model over the grid-aligned method for 
second-order computations of subsonic airfoil flows such as this. 


8.2 Navier-Stokes Computations 

Only second-order spatial accuracy is used for all Navier-Stokes computations. 
No limiting of higher-order terms is employed. All solutions are advanced in time 
using the implicit approximately-factored algorithm, and convergence is assumed 
to be reached when the L 2 -norm of the residual drops by at least 4 orders of 
magnitude and/or the lift and drag values settle down and do not vary significantly 
with further iterations. 


8.2.1 Shock /Boundarv-Laver Interaction 


A shock/boundary-layer interaction is computed over a flat plate in a domain 
that is 1.84 units wide by approximately 1.29 units high. The plate is assumed 
to start 0.24 units from the left. The finest grid employed is a 93 x 141 with 
minimum spacing at the wall of 9 x 10~ 5 units. It is shown in figure 8.49. Coarser 
meshes of 47 x 71 and 24 x 36 are obtained by removing every other point from 
the next finest mesh. 

The inflow Mach number is 2.0, and the Reynolds number per unit length is 
296,000. Laminar flow is assumed. A shock wave enters from the left of the domain 
at a height of 0.7930 units. It has an angle such that it impinges on the plate 
1.0 unit from the leading edge. The flow is turned through an angle of 3.1°. The 
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Figure 8.49: Shock/Boundary-Layer Interaction 93 x 141 Grid 
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boundary conditions at inflow (below the shock) are: p = 1, pu = 2, pv = 0, pE = 
3.78571. Above the shock and along the top wall: p = 1.13038, pu = 2.18509, 
pv = -0.11834, pE = 4.23855. At the right wall all variables are extrapolated 
from the interior using second-order interpolation. Along the bottom wall in front 
of the plate symmetry conditions are applied, and on the plate itself no slip, zero 
pressure gradient, and adiabatic wall boundary conditions are assumed. For this 
problem, the thin-layer Navier-Stokes equations are employed: it is assumed that 
viscous terms arising from derivatives in the streamwise direction are negligible in 
comparison to the viscous terms arising from derivatives normal to the wall. 

Pressure contours for the grid-aligned model on the three grids are shown 
in figure 8.50, while 5-wave model results are shown in figure 8.51. (The 5-wave 
model results are obtained by restarting the grid-aligned solutions with 9 d frozen.) 
Values in the figures are plotted from 0.84 to 1.44 in steps of 0.03. Results are very 
similar for both methods, with the exception that the 5- wave model gives slightly 
sharper shock resolution. Also, some oscillatory behavior is evident near the plate 
leading edge as well as at the shock inflow location for the 5-wave model. Skin 
friction plots over the plate surface are given in figure 8.52 for the two methods. 
In figure 8.52(a) it is seen that the grid-aligned method gives reasonable levels 
on all three grids, with the greatest discrepancies between the solutions occurmg 
in the separated region and near reattachment. The 5-wave model, on the other 
hand, gives very high skin friction levels on the coarsest mesh (as shown in figure 
8.52(b)), indicating that the velocity profiles are predicted to be too full. This can 
be seen in the computed velocity profiles at two locations on the plate in figure 
8.53, in comparison with grid-aligned results on the same grid. On the two finer 
grids, the 5-wave model gives consistent skin friction results in good agreement 
with the grid-aligned model. 

It is not known why the 5-wave model predicts fuller profiles on the coarsest 
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mesh for this case. Aside from this discrepancy, the model behaves much as 
expected, considering the results from the earlier Euler shock reflection case. For 
this second-order computation, the shock structure is captured slightly better 
than when the grid-aligned scheme is used, but some oscillations also appear in 
the flowfield. The very small increase in shock resolution does not seem to be 
worth the added expense and increased oscillatory behavior of the 5-wave model. 


8.2.2 Subsonic Separated Flow Over an Airfoil 


A specific case where an advantage of the 5-wave model over the grid-aligned 
model is fully realized in a second-order computation is for viscous separated flow 
over a NACA 0012 airfoil at M = 0.5, a = 3°, and Re = 5000. Full Navier-Stokes 
computations are performed on a series of C-meshes. The finest, a 257 x 97 grid 
with outer boundary extent of 14 chords and average minimum spacing on the 
body of 2 x 10 -4 chords, is shown in figure 8.54. There are 176 cell faces on the 
airfoil. Coarser meshes are obtained by removing every other gridpoint from the 
next finest mesh. Subsonic inflow /outflow boundary conditions are applied on the 
outer boundary, and no slip, zero pressure gradient, and adiabatic wall conditions 
are used at the body. Continuation conditions are applied on the wake cut. 

This test case was first discussed by Venkatakrishnan [31]. He found that 
on even reasonably fine meshes, the grid-aligned upwind scheme does a poor job 
for this case since there is a detached shear layer emanating from about mid- 
chord on the airfoil upper surface which is not oriented with the grid. The shear 
is misinterpreted by the grid-aligned model as a combination of shear and com- 
pression/expansion, with the end result of a distortion in the computed pressure. 
Grid-aligned results on the three grids showing this pressure distortion are shown 
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Figure 8.54: NACA 0012 257 x 97 C-Grid 
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in figure 8.55. Pressure contours axe plotted in increments of 0.005. As the grid is 
refined the distortion diminishes in extent, but does not goes away on these grids. 

The 5-wave model is run by restarting the grid-aligned solutions with 6 j 
frozen. Results are plotted in figure 8.56. On all three grids, the pressure distortion 
is essentially eliminated. Lift coefficient values, plotted as a function of the inverse 
of the number of gridpoints in figure 8.57, show that both models are approaching 
the same value of c/ « 0.055 as the mesh is refined. In spite of the fact that the 5- 
wave model does a good job of reducing the pressure distortion on all three meshes, 
the lift is seriously overpredicted (by about 120%) on the coarsest one. This can 
be seen in a plot of surface pressure coefficient on the 65 x 25 grid, given in figure 
8.58. The surface pressure values given by the 5-wave model are significantly lower 
than the grid-aligned values over the upper surface, as well as slightly higher over 
the lower surface. On the other grids the two models are in closer agreement, 
although the grid-aligned model is always closer to the “correct’ extrapolated lift 
coefficient. Drag coefficient values are plotted in figure 8.59. Figure 8.59(a) shows 
the total drag coefficient, while figure 8.59(b) breaks the drag coefficient into its 
friction drag and pressure drag components. Both models do reasonable jobs of 
predicting the drag on the two finest meshes. 

In summary, the 5-wave model can essentially eliminate pressure distortions 
which arise in a separated viscous flow computation due to misinterpretation of 
oblique shear waves by the grid-aligned scheme. However, the surface pressure 
levels (and hence the lift coefficient) predicted by the 5-wave model are more 
sensitive to grid density than those predicted by the grid-aligned model, so a 
fairly fine mesh is required for a reasonable level of accuracy in this regard. Since 
the grid-aligned method shows pressure distortions on even the 257 x 97 grid, the 
5-wave model appears to be worth the additional effort in this case. 

All results obtained for the Navier-Stokes cases to this point using the 5-wave 
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Figure 8.57: Lift-Coefficient Grid-Convergence Study for Viscous Airfoil- 

Flow 
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Figure 8.58: Viscous Airfoil-Flow, Surface Pressure Coefficients, 65 x 25 
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a) Total Drag Coefficient 
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b) Friction and Pressure Drag Coefficients 
Figure 8.59: Drag- Coefficient Grid-Convergence Study for Viscous 


Airfoil- Flow 
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model have employed ^-limiting. When 6 ' d -limiting is not employed, the odd-even 
point decoupling mentioned in Chapter 7 can occur in the boundary layer. An 
example is shown in figure 8.60. Pressure contours over the rear half of the airfoil 
are shown for computations on a 129 x 49 grid. Several of the contours show the 

oscillatory behavior. 



CHAPTER 9 


EXTENSION TO THREE DIMENSIONS 


9.1 Governing Equations 


The nondimensional Navier-Stokes equations in curvilinear 
conservation form can be written in three dimensions as: 


where 


au* dF* dG* dH* dF* dG* v dH 

~dT + ~dl + d v + 8C d£ + drj + dC 
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F* = - 

*• V T 


G* = - 


h: = - 


0 

£s r ll + £y T 12 + tz T 13 

i x T 2\ + iy T 22 + £z r 23 

| £xT 3 i + £y T 32 + £z r 33 
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0 

Cr T ll + Cy r 12 + Cz r 13 
Cx t 21 + Cy T 22 + Cs r 23 
Cx t 31 + Cy r 32 + Cz r 33 

CJViT 1 i - Ql) + (y{VjT 2 j - Qi) + Cz{VjT 3 j - Q3) 


.Vz 


(9.6) 


(9.7) 


(9.8) 


u * = txU + tyV + £* w 

v* = Vx u + T] y v + TJzW ( 9 - 9 ) 

W* — Cz u + Cy u + Cz^- 

All variables are nondimensional as defined in Chapter 2. The terms u , v , and 
w* are the contravariant velocity components, and V u V 2 , and V 3 represent u, v, 
and in, respectively. T i; and Qi are given by (2.14) and (2.15), but now 


0 

dX 1 

8 

dX 2 

d 

dX 3 


d d d_ 

= t*dt +Vx dr, +Qx d<; 

. d d . d_ 

-t* dl +T,y d v +Cy d( : 
a d . d_ 

~t*di +r]z dr ) +Cz ac* 


( 9 . 10 ) 


The ideal gas equation of state is given by 


( u 2 + v 2 + in 2 \ 

V = (7 - l )P [E ~ 2 ) ' 


(9.11) 


Also, Stokes’ hypothesis, A + (2/3)/* = 0, and Sutherland’s law (2.34) are assumed 
to apply, and 7 is taken as 1.4 and Pr is taken as 0.72. 
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In three dimensions, the inverse of the metric Jacobian J is given by: 

J -1 = X{(y v z( - y^z v ) - y^(x v z <: - x^z v ) + Z{(x v y <; - x^y v ), (9-12) 

and the metrics of the transformation are 

£. = J(yr,z c - y f *„) (y = J(2r,x c - z c x v ) ( z = J{x v y< - x f y„) 

Vx = J{y<; z t - Vi z <) Vy = J( z < x t ~ hh) Vz = J{ x <;yt - x ty<i) (9.13) 

Cx = J(yt 2 v ~Vv z t) C v = J ( z t x v - z v x t) (z = J( x ty v ~ x vyt)- 

The Euler equations are the same as the Navier-Stokes equations (9.1), except 
that the viscous terms F* , G* , and H* are taken as zero. 


9.2 Spatial and Temporal Discretization 


Since only Euler computations are performed in three dimensions in this 
study, the discretizations described here will not include viscous terms. The basic 
scheme employed is CFL3D [32]. It is an implicit finite-volume method, in which 
the left-hand side is approximately factored and A t U^ n ^ is solved for in three 
sweeps through the mesh: 
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+ *« 
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+ *< 
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dV 

dG* 

dU 
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dll 


A t U' = T 


A - u " = JKt A,v ' 


A * u< “ ) - jAt ‘ A<u ”’ 


(9.14) 


where the terms A t U' and A t U" are intermediate results. The conserved variables 
are updated at the cell centers using 


u (n " t ' 1> = U (n) + A,U (n ). 


(9.15) 
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The right-hand side T in (9.14) is defined by 

T = -(« £ F* +f,G* +«<H'), 


(9.16) 


and 




(9.17) 


The right-hand side terms in (9.17) are the fluxes at the six faces of a computa- 
tional cell, evaluated through the use of a dux function, as described in sections 
9.3 and 9.4. The right-hand side T can also be written in the form 


T = - y ‘lyA.l/, 


(9.18) 


1=1 


where the summation is over the six faces of the computational cell, Ax4* is the 
area of cell face t, and 3> is the normal flux per unit face-area, given by: 


$ = 


Pig 

pq g u + p{c x ) g 
pq g v + p(c y ) g 
pq g w+p(c z ) g 

pq g H 


(9.19) 


The variables c„ c y , and c 2 represent the components of the unit-normal direction 
vector n in the x, y , and z directions, respectively. They are written in spherical 

coordinates as: 

C x = COS0COS^> 

Cy = sinflcos^ (9.20) 

c z = sim/>. 

The subscript g is used to indicate the grid-face normal angle. For example, 
(c) p = cosOg cosip g . The angles 6 g , g and the grid-face normal direction vector 
n g are pictured in figure 9.1. The variable q a is the outward velocity normal to 

the cell face, given by 


q g = u(c x ) g + t>(cj,) fl + w(c z ) g . 


( 9 . 21 ) 
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The implicit spatial derivatives of the convective and pressure terms are spa- 
tially first-order accurate, resulting in block-tridiagonal inversions for each sweep. 
Approximate left-hand side Jacobians of the form given in Chapter 3 are employed 
regardless of whether the grid-aligned or the grid-independent model is employed 
on the right-hand side. 


9.3 Grid- Aligned Flux Function 


The normal flux per unit face-area in three dimensions can be computed using 
the grid-aligned flux function: 


4 

# = \ ($l + $r) - ~ y^|A fc |fl fc Rfc, 

^ 1 *:=i 


where the waves are represented by the vectors 
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(9.22) 


(9.23) 


(9.24) 


(9.25) 


R.4 = 


1 

u 

V 

U) 


L \{u 2 +V 2 + W 2 ). 


(9.26) 
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The terms Ar x , A r 2 , and Ar 3 are given by 

An = {(c x ) 2 - 1) Au + ( Cl )(c y )Av + (c x )(c z )Ain 

Ar 2 = (c s )(c z )Au + ((c y ) 2 - 1) Au + (c s )(c 2 )Auj (9.27) 

Ar 3 = (c z )(c z )Au + (c z )(c y )Au + ((c z ) 2 - 1) Atu 

and the subscript g indicates that the components ((<:*)„, (c y ) g , (c z ) g ) of the grid- 
normal direction n g are used. The Roe-averaged outward velocity normal to the 

cell face is 

q g =u{c x ) g +v{c y ) g +w{c z ) g . (9-28) 

The vectors (9.23) through (9.26) represent, respectively, + acoustic, - acoustic, 
shear, and entropy waves. R x , R 2) and R 4 are eigenvectors of the linearized three- 
dimensional Euler equations, while R 3 is a linear combination of the remaining 
two eigenvectors, which, along with R x , R 2 , and R 4 , span the eigenspace of this 
system. The expression for the combined shear wave R 3 is derived in Appendix 
B, and a geometric interpretation of its effects in state space is given below. 

The vector of wavestrengths is given by: 

^(Ap + pa Aq a Y 

^j{Ap-paAq g ) (9.29) 

£ ’ 

i (a 2 A°p - a p) . 

where 

A q g = A u(c x ) g 4- A u(c„) 9 + Aw(c z ) g . (9-30) 

The wavespeeds are all in the n ff -direction and are given by: 

\i=q g + a 

^2 = Qg ~~ ® 

A3 = q g 

A 4 = 9s* 



(9.31) 
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The Roe-averaged (hatted) values are still defined by (4.9). 

In three dimensions, the (Au, Au, Ap)-space interpretation of the wave de- 
composition is not as helpful as it was in two dimensions, since the jump Au; is now 
important. Instead, when describing the action of the three-dimensional waves, 
the Ap- “direction” is ignored, and jumps are plotted in (Au, Au, Au;)-space. 
Hence the change in pressure caused by an acoustic wave cannot be represented. 

Figure 9.2(a) shows the change in velocity associated with a +n g acoustic 
wave, while figure 9.2(b) shows the change for a -n g acoustic wave. Each of these 
waves also causes a positive change in pressure and a change in density, which are 
not pictured. Similar to the two-dimensional case, the entropy wave causes only a 
change in density (with no change in velocity or pressure), so it is not representable 
at all in (Au, Au, Aiu)-space. Finally, there are two independent eigenvectors, 
as discussed in Appendix B, which represent shear in the plane perpendicular to 
the n g -direction. Any two eigenvectors in this plane are independent provided 
that the directions of the corresponding velocity jumps are perpendicular. An 
example of two such waves is given in figure 9.3. These waves propagate in the 
redirection, which is normal to the change in velocity across both of the waves. 

The two shear waves can be replaced by the single shear wave (9.25) with 
strength p/a. Its net result is the sum of the velocity changes across the other 
two. This single shear is, in effect, the wave that causes a change m states from 
R to R\ where R' is the point in (Au, Au, Atu)-space on the fine from L with 
direction n g where the perpendicular line from R intersects. This is illustrated in 
figure 9.4, along with hypothetical + and — acoustic waves. 
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9.4 Grid-Independent Flux Function 


The dominant direction of wave propagation for the grid-independent wave 
model is again assumed to be the velocity-difference direction. For three-dimen- 
sional flow, this direction is given by defined by the angles (Och^d), which are 
obtained using: 


0d = tan 


rpd — tan 


-(=) 

_j / Aw * sig^Atf) ^ 

\ \fKu* + Au 2 / 


(9.32) 


The angles are each defined from — ^ to 

As in the two-dimensional grid-independent model, the wave types and stren- 
gths for the three-dimensional model describe the difference between the left and 
right states. The wave-propagation direction can be frozen as n' d (defined by the 
angles (^,V»d))» resulting in the requirement of a shear wave of type similar to 
the shear wave used in the grid-aligned model. The end result is, again, a 5-wave 
model, which is computed using: 


$ = 


-(*L + *a) 


jEMa***- 


(9.33) 


The Rjt are given by 


Ri 


1 

u -I- a(c x )' d 
v + a(c y )' d 
w + a(c z )' d 
H + aq d _ 


1 


u - a{c x y d 
v - a{c y y d 
w - a(c z )' d 
. H-aq' d . 


(9.34) 


R 2 = 


(9.35) 
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R-4 — 


Rs = 


where 


0 

—a( c x)d 

R3 = I ^{ c y)d 
-a( c *)'d 
-aq' d J 

1 

u 

v 

w 

L|(u 2 + v 2 +W 2 ) 

0 

—a(Ari)' d 

-o( Ar 2 )' d 

-a{Ar 3 )' d 

— d{u(An)d + v(Ar2)J t + ^(An)^} 


q d = u(c x )' d + v{c y )'d + u>{cz)'d- 


(9.36) 


(9.37) 


(9.38) 


(9.39) 


The parameters c„ c y , c„ and An, Ar 2 , Ar 3 are still defined by (9.20) and 
(9.27). The subscript d and the prime (') indicate the frozen direction n' d . For 

example, ( c x )' d = cosfljcosi/^. 

The vector of wavestrengths is given by: 


n = 


rfg + ^(Au{c 9 y d + Av( Cy y d + Aw{c z y d ) 
- f| (A u{c x y d + At?(c y ) ^ + Ain(c z )J i ) 
(/? — 1)| (Au(c*)Jj + Av(c y )' d + Aw(c z ) d ) 
£ (a 2 A p - A p) 


R 

a 


(9.40) 


The parameter (3 is similar in form to that derived in two dimensions. For the 
minimum-pathlength model, 

Ap/{pa) 


/ 3 = min 


,1 


(9.41) 


I Au(c a: )J i + A’u(cj,)' (1 + Aw(c z )' d 
The minimum-area model, used for all computations in this study, results when 

Ap/ {pa) . 


mm 




A u(c x )' d + Au(c v )' (i + Aw(c z )' d 


!l - ■ 


(9.42) 
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The parameter (3 is generally limited to be no less than 0.05 and is frozen along 
with n' d as an aid to convergence. 

The effects of the acoustic waves Ri and R 2 as well as the shear wave R 5 
can be plotted in ( An, Au, Aio)-space. Plots look the same as the grid-aligned 
wave plots (figures 9.2(a), 9.2(b), and 9.4) except that the angles are now given 
by 0' d and The shear wave R 3 is henceforward referred to as the {n' d + f) 
shear wave, in order to parallel the {d' d + f ) shear wave from the two-dimensional 
theory of Chapter 5. It causes a change of velocity in the n^-direction, as pictured 
in figure 9.5. Therefore the direction of propagation of the wave itself must be 
90° to this. Unlike the two-dimensional case, however, in three dimensions there 
are an infinite number of directions perpendicular to the n^-direction. This is 
tantamount to saying that the difference in velocity between two states L and R 
(with no pressure difference) can be attributed to an infinite number of shear wave 
types, each of which travels in a different direction perpendicular to n' d . 

An illustration of this concept is given in figures 9.6(a) through (e). For 
this example a first-order computation is assumed, where two neighboring cells 
have velocities V L = (1,0,0), F* = (1,1,0). Also, the ^-direction is assumed to 
be computed from these velocity values (i.e. they are not frozen from an earlier 
time step). Figure 9.6(a) shows the geometry, including the velocity-difference 
vector AV which indicates the direction n d . One type of shear wave that could 
describe this difference in states is a “layer-type” shear, as pictured in figure 9.6(b). 
Across this shear plane the velocity vector both rotates and lengthens. There is 
no velocity component normal to the “layer-type” shear wave, so it is steady with 

drift velocity us = 0. 

Other interpretations of the velocity difference are shown in figures 9.6(c) 
through (e). In each of these, the shear wave does have a drift velocity u s . This 
drift velocity is always normal to the redirection. It is a maximum in figure 
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Figure 9.5: Effect of (n' d + f ) Shear Wave 
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9.6(e), for the case of the “2-D-type” shear wave. It is termed “2-D” because the 

— ♦ — ♦ 

drift velocity is parallel to the plane defined by the velocity vectors V L and V R 
(with origins of the vectors at the same point). Across a “2-D-type” shear wave 
the velocity vector lengthens, but there is no rotation. It should also be noted 
that although the wave passes through both points L and R in figure 9.6(e), it is 
assumed that L lies infinitesimally to one side of the wave and R to the other. 

The initial strategy used to determine the direction of propagation of the 
(n' d + |) shear wave involved the grid geometry. The distance vector D g , which 
connects the cell centers of the neighboring states L and R (see figure 9.7), would 
be compared to the plane defined by the velocity vectors. If D g was at right 
angles to the velocity plane, then a “layer-type” shear would be assumed. This 
corresponds exactly to the situation pictured in figure 9.6(b). If D g lay parallel 
to the velocity plane, then a “2-D-type” shear would be assumed. Anything 
else would fall somewhere between the two types of shear. An easy method of 
implementation would be to find the normalized projection of D g onto the plane 
which is perpendicular to the direction vector n' d . This direction would then be 
taken as the direction of shear wave propagation. The shear wave drift velocity 
would simply be the component of (it, ■£? , ti>) taken in this direction. 

However, this strategy has two undesirable properties. First of all, it brings 
grid-dependency back into the solution. Second, and most importantly, it does not 
allow for sharp capturing of steady three-dimensional oblique shear waves through 
which the velocity undergoes rotation. For example, the steady-state type of 
shear shown in figure 9.7 is not interpreted by this model as pictured; instead, the 
difference in states is construed to be due to a different shear wave with non-zero 
drift velocity. This misinterpretation adds dissipation to the computed solution 
and smears the final steady-state result. 

The best propagation direction for the ( n' d + f ) shear wave seems to be the 
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one that yields the lowest dissipation. This is given by the “layer-type” shear 
wave interpretation, which always has drift velocity us = 0 when n<* is computed 
every time step. When the direction n d is frozen, the drift velocity is no longer 
necessarily zero, but it is generally very small. The unit normal direction vector 
is computed as the direction perpendicular to both Vj, and Vr using the following 
formulas: 

a = vlwr — wlvr 
b = wlUr — u L w R 

(9.43) 

C = ULVR — VLUR 

d = y/a? + 6 2 + c 2 . 

Then, if d ^ 0, 

((«.)».(<i)».(«*)s)= (3.3.5)- < 9 ' 44 > 

If d = 0 the arbitrary unit-normal direction (which is perpendicular to the re- 
direction) is chosen: 

((c«)s,(c v )s,(c*)s) = (-cos9 d sinil> d ,-sin0 d sinj> d ,cosij> d ) . (9-45) 

The drift velocity of the shear is then computed as the component of (u,v t w) in 
the ((c.) 5> {c y )s, (c z )s)-dmection: 

U S = fi(c x )s + v(c y )s + tt)(c z )s. (9.46) 

Again, if the n^-direction is not frozen, then ((c x )s, (c y )s> (c z )s) is computed 
from the same V L , Vr that formed ( u,v,w ), and u s turns out to be identi- 
cally zero. However, if the n^-direction is frozen from a previous iteration, then 
((c*)s,(c y )s,(c*)s) is frozen as ((c x )' 5 ,(c y )' 5 ,(c z )' s ) and u' s = u(c x )' s +u(c y )' 5 + 
ttj(c z )' 5 may have a finite value. 

Finally, the components of the wavespeeds of all five waves in the grid-normal 


direction are written: 


Ai = ( q ^ + Q-){( c x)d( c x)p + ( c y)d( c y)y + ( c *)d( c *)y} 

A 2 = (?a “ fl){( c *)d( c *)y + ( c y)d( c s/)y + ( c *)d( c *)y} 

A 3 = n' s {(c x )' 5 (c x ) 9 + (c y )'s(c y ) 9 + (c z )' 5 (c 2 ) 9 } ( 9 - 47 ) 

A 4 = <ld{{ c x)'d{ c x)g + i c y)d( c y)g ( c *)d( c *)fl} 

A 5 = q'd{{ c z)'d{ c *)fl + ( c y)d( c v)$ + ( c *)d( c *)ff}- 


9.5 Stability Analysis 

Following similar linearization procedures as for two dimensions, as well as 
making the assumption that all cells are cubes with edge length As, the three- 
dimensional Euler equations can be written as 

Ai ®U = -LU, (9.48) 

at 

where L is a matrix of difference operators. For first-order spatial differencing, 
the Fourier transform of the right-hand side turns out to be: 




— cos£^)+ 

- cosC (fc) ) + 


(9.49) 
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For second-order differencing it is: 

5(-L) = jw ~JkT~ [|^ > (”s ) )| (- 2 cos C (l) + 

{a >(n g ' ) + u(rig ) + u}(n g )} L 

icos2C« + |) +i(|§(4°)) ( 2si "C <i) - r in2C<i, ) + 

\t^)\ (-2cosC 0) + ^cos2C (i) + 

I) + i(^(4 ;) )) ( 2 “< (i> - 5 ri “ 2 < (<> ) + 

|d(4‘>)| (-2cos<<‘> + icos2C (fe) + 

l) +i (^( ,i< *‘ ,) )( 2sin <“ ) -r in2<( ‘ ) )]' 

The variable u is the CFL number, defined as 

v = {w(n ( 9 l) ) + u;(n^) 


(9.50) 


(9.51) 


for cubic cells with edge length A s. The rig \ n^\ and ^ are the direction- 
vectors normal to the grid faces in the i, j , and k directions, respectively. They 
are defined by the angles ff g and ip g (see figure 9.1). In the present analysis, 6 g and 
are prescribed for the direction Tig \ then iig^ is taken as ( 0 g ^pg d" 7r /2). The 
vector n 9 ^ is taken as the vector perpendicular to the other two. Also, w(n 9 ), 
u>(n^), and u ;(n 9 ^) are the maximum wave speeds |g| + a in each of the grid 
directions. $( — L) is a complex- valued 5x5 matrix. 


The three-dimensional flux Jacobian matrix is given by d&/dU — 

0 (c*) 9 (c 9 ) 9 ^ (c z ) g __ 0 

(Cx)g4>-'uqg {Cz)gbu + q g ( c y)g^-( c x)gf v ( c z)g u ~( c x)gf \ c x)gf 

\c y )g(j)-vqg ( C x )gV-(Cy)gfU (Cy)gbV + q g ( C Z )g V ~ ( C y ) g /' W (Cy) g f 

( C z )g<p — Wq g ( C x )gW — (C z )gfu ( C j, ) g W — ( C z ) g fV ( C Z ) g bu> + ^g { C z) g f 

(2<^ — 7 E)q g Cm C43 ^44 79s . 

(9.52) 


where, for brevity, 

^42 = (cx)p(7^ - <t>) ~ fuq g 

C43 ={ c y ) g{-fE - <i>) - fvq g 

C44 = (c z ) 9 (7-E fwq g . 


(9.53) 
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Also: 


*3(2-7) 

/ = (t - 1) 

<j> = (u 2 + v 2 + w 2 ) , 


(9.54) 


and g 9 is defined by (9.28). 

|D| is the matrix that satisfies 


D 


AU = HfcR-fc* 


(9.55) 


For the 5-wave model in three dimensions, the elements of |D| are given by: 


D 

D 

D 

i 

Id 


j i 


J'2 


|A 4 (R-4)j + <£6 - 6* - - 6™ + 

A5 {(pi)j« + (p 2 )iU + (P3)» 

= - (7 - i)u6 + 6 - 1^5 1 (pi )i 


= - (7 - l)i)6 + 6 - 1-^5 1 (P2)i 
j3 

= - (7 - l)™£l + 6 - 1^5 1 (P3)i 

j4 


jS 


=(7 - 1 )6, 


where 


(R 2 ) i -2|A,|(R 4 ) i } 

(ROi) 


<• = 2 ?{M<^ + 

6=^(c,)i{|^i|(Ri)i 
«4=^(c=)i{J^J(R0j-|^|(R,)i} 


and 


(pjfe)i = 0 k = 1,2,3 

(pi ) 2 = ( c x)d — 1 (P2 )2 = { c x)'d( c v)d 

(Pl)3 = (Cj ,)'d( C x)'d Ms = (Cy)d - 1 

(Pl)4 = (c JB )rf(c a5 )J i (p 2 )4 = ( c z)d( C v)d 

(pjfe) 5 = (Pl)*+1^ + Mk+lV + (p S )fc+lt& 


(Ps )2 = ( c x)'d{ c z)'d 
(Ps)s = i c y)di C z)d 
(P3)4 = (c.)2 ~ 1 
k = 1,2,3. 


(9.56) 


(9.57) 


(9.58) 
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In obtaining these elements of the |D| matrix, it is assumed that the wavespeed 
of the ( n' d + -|) shear wave is zero. 

Although no three-dimensional computations using explicit time-marching 
were done in this paper (only the modified implicit code CFL3D was used), a 
linearized stability analysis for the explicit scheme is included here for complete- 
ness. The analysis of the 4-stage time-marching scheme (6.18) is done in the same 
way in three dimensions as in two. However, since there are now three perturba- 
tion wavelengths < (i) , and C (fe) , each is cycled through 9 values from -tt to 
7 r inclusive, for a total of 729 conditions. Five eigenvalues are obtained at each 
condition. Since there are now more independent parameters to vary (the flow 
direction ny, the grid-normal direction n s , and the wave-propagation direction 
n' d each are defined by two angles 9 and if}), it is even more difficult to perform 
a thorough numerical analysis. Through extensive numerical experimentation, a 
“worst-case” Fourier footprint (of all variations tested) has been determined for 
first-order spatial differencing with the parameters M=100, 6/ = 0°, iff/ = 0°, 
(3 = 0, = 0°, = 0°, 9' d = 22.5°, and if>' d = 22.5°. 

Plotting this “worst-case” Fourier footprint along with the stability boundary 
of (6.18) with r] = 0.15 in figure 9.8, it is seen that the maximum u that yields sta- 
bility of the three-dimensioncil explicit scheme is about 1.55. This is slightly lower 
than the stability limit for two dimensions. For second-order spatial accuracy, the 
maximum v is about 0.77, as shown in figure 9.9. It is expected that for lower 
Mach number flow these CFL number restrictions would be relaxed somewhat. 

For three-dimensional implicit time-marching, the generalized eigenvalue pro- 
blem 

[3(M) + 3(-L)]x = g[S(M)]x, (9.59) 

arises, where g is any of the complex eigenvalues, and [Sj(M) + S(— L)] and [^r(M)] 
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are cor 


jmplex 5x5 matrices. For the three-factor approximately-factored left-hand 
side, (9.14), the Fourier transform of M is given by: 


3(M) ={ I + 

dl 
.a U' 

1 + 

.at 

i + 



>) + w (4") + 
_ \ 

^(4 fc) )} L 

\&U 

:(4°)) 

|sin£^ 

)• 





V 


a# 
1 

M4' 

+ u; ( n p^) + 


dU 


j sin£^ 






V 

I 

’i d& 


(9.60) 


{u;(n^) + 01(71^) + u>(nj ; )} 


|J(n< fc) )|(l- co sC (fc) )+ 


when the approximate left-hand side Jacobians of the form given in Chapter 3 are 


employed. 

The stability characteristics are determined by cycling through 9 of each of the 
frequencies £ (i) , C (j \ and £ (fc) from 0 to 2?r. For first-order spatial differencing, the 
same “worst-case” parameters as those used for the explicit analysis were found to 
give the strictest stability limit out of the many variations tried. These parameters 

are: M=100, 9j = 0°, = 0°, tf = 0°, # = 0°, 9' d = 22.5°, and ft = 22.5°. 

Variations in /3 were found to make very little difference at these conditions. A 
plot of the maximum eigenvalue, average eigenvalue, and smoothing factor is given 
in figure 9.10. The maximum CFL number allowable turns out to be about 1.5. 

For second-order spatial differencing, the same “worst-case” parameters were 
again used, but this time (3 was found to have an effect on the results. When 
(3 — 0, the analysis shows the scheme to be unstable except below v = 0.005. 
This is shown in figure 9.11(a). When (3 = 0.05, figure 9.11(b), the stability limit 


increases to about 0.3. 
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Figure 9.10: Eigenvalues as Function of CFL Number for Implicit Scheme 
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Based on the results from two dimensions, it is believed that these stability 
restrictions for first and second-order can be relaxed somewhat for lower Mach 
number flows. It is also likely that the stability properties of the 5- wave model in 
three dimensions can be improved dramatically by modifying the left-hand side 
approximate Jacobians as described at the end of Chapter 6, as well as solving 
bio ck-p e nt adi agon al systems for second-order computations. These modifications 
were not attempted in the current study. 


9.6 Monotonicity Analysis 


The monotonicity analysis for the three-dimensional Euler equations proceeds 
much the same way as for the two-dimensional equations. The monotonicity 
property is insured when 


e.v. 


(' an +U ~ ) ~° e - v - - °’ (9 - 61) 

V^Ui+i j t ii J \^Ut_i ,j,k J 


where e.v.(-) represents “the eigenvalues of (•)”, and the grid-normal n g , defined 
by the angles (Og,i^g), is varied over the full range of possible angles. The 5-wave 
model can again be written in a form similar to (7.9), but it is assumed for this 
analysis that the wavespeed associated with the (n' d + f ) shear wave is identically 
zero (which would be exactly true if the wave-propagation direction was never 
frozen). Therefore this shear wave is ignored in the analysis. All remaining waves 
travel in the n^-direction, so the flux function per unit face-area can be written: 

+ *«,.) - |[R*]|[A*]cos(ni - n,)|n*. (9.62) 


There are only four waves which are assumed to travel in the ^-direction, given 
by (9.34), (9.35), (9.37), and (9.38). However, in order for all the matrices in 
(9.62) to come out as 5 x 5 matrices, it is necessary to employ five waves in R*. 
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Therefore the combined shear wave R 5 of (9.38) is broken back into two shear 
waves, chosen arbitrarily, which are perpendicular in ( Ati, Au, Aw)-space. 


R 5 a — 


0 

— asinflj 

acos8' d 

0 

a(— usind' d + ucosflj) J 


(9.63) 


Rsi, — 


0 

— acos9' d simp' d 
— asin^sinV»j 
acosV'Jj 
af!, 


(9.64) 


where 


r' d = —ucos8' d s\nifi d — vsinO d siml)' d + wcosijj d . 

The wavestrengths associated with each of these waves are: 

fi 5o = t ( — Ausin^ + Ancosflj) 

a 

fl 5b = ^(-Aucos^sin ^ d - Avsin^sin^ + Awcos ^ d ). 


(9.65) 


(9.66) 


Hence the [R*] matrix is made up of R.1,2,4 from (9.34), (9.35), and (9.37), as 
well as R 5a and R 56 . The matrix of wavespeeds in the redirection is [A*] = 
diag (q' d + a, q' d - a,q' d ,q' d ,q' d ), where q' d is defined by (9.39). Making the same 
assumptions as was done in Chapter 7 , the linearized analysis gives 


dU i+1>i , fc 2 8U 


1 S * («») - i[B.'l|[A']cos(iSi - n„)| a0 


8 U r 

18 *,,, 1 , 6 . 11 . J.i ^ _ * ll^Ol 
au,-,.,., = 2 8 U (ng) ~ 2 [R1||A |cos( ' 9,1 au t ' 


(9.67) 

(9.68) 


The derivative matrix in (9.67) is given by d£l*/dXJn. — 

<P f-u 1 P(cx)'i fv | At _|_ 

2d* ~ it ~2d* + 2 a 2d* + 2a 2d* + 2 a 

_ „ i - . si/ \ t r * «[/•_!. 


££. 

2a 

1-4 


JL + 

2a 5 T 


4 U . 

/u HsxYj, ix _ PLsih. 

2d 3 2a 2d 7 2a 2d 5 2a 

IX II L^S. 


a 

4 

d 


sin#!, 


d 3 

cosflj 


u 

cos^sin- 0 j 


a 


0 

COS- 0 J 


2a 3 

2d 3 

_JL 

d 2 

0 

0 


(9.69) 
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where the variables / and <f> are defined in (9.54), and: 

t' d = — usinf^ + vcos0' d (9.70) 

The derivative matrix in (9.68) is given by (<?n*/t?Ui,) = — (5rt /5 Ur). 

For the three-dimensional monotonicity analysis, the Mach number M, flow 
angles and /3 are chosen, then 9 g , i> g , 9 d , and ip' d are each varied inde- 

pendently from —90° to 90°, with incremental change of 7 t/ 8. Eigenvalues are 
computed for each condition. If they meet the criteria of equation (9.61), then 
monotonicity is preserved at that condition. The results are plotted as allowable 
| n' d - n { \ vs. \n g — n/|, for various <£„, where <j> p is the angle between the plane 
defined by the vectors (n g ,nj) and the plane defined by the vectors {n d ,nf). The 
quantities within the absolute value signs indicate an angular difference between 
two vectors. When plotted this way, results are independent of the flow vector n f . 

A result is given in figures 9.12(a) through (c) for M = 3, and /3 = 1. Figure 
9.12(a) shows results for the case when <f> p < 15°. The allowable region for mono- 
tonicity looks very similar to the two-dimensional region at these same conditions 
(see figure 7.2(a)). It includes the grid-aligned wave model n' d = n g . Notice 
that for the three-dimensional case, in contrast to two dimensions, the absolute 
value of the angular differences are plotted so that only positive differences are 
given. This is done because of the difficulty associated with assigning a positive 
or negative angular difference in three dimensions. When 15° < <f> p < 30°, the 
plot of figure 9.12(b) results. Here, the monotone region is similar to that in fig- 
ure 9.12(a) except that the grid-aligned model is no longer representable. When 
30° < 4> p < 45°, the monotonicity region diminishes significantly in size, as shown 
in figure 9.12(c). Finally, when 45° < <f> p < 90°, then no region is monotone, 
according to this analysis. 

A specific example is taken from this case. Referring to figures 9.12(a) through 
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(c), when \n g — n/| « 75°, the allowable \n' d -n f \ goes from about 25° to 75° for 
(j) p < 15°, from about 25° to 65° for 15° < <f> p < 30°, and from about 35° to 45° 
for 30° < 4> p < 45°. A sketch is first drawn in figure 9.13(a) of the nj vector and 
the ng vector, separated by 75°, with the allowable n d directions in the (n.g,7iy) 
plane (cf> p = 0°) indicated by shading. Next, in figure 9.13(b), the allowable n d 
directions in all three dimensions are indicated by including the results from the 
cases when the (n^,ny)-plane differs significantly from the (n s ,ny)-plane. 

A second case using M = 0. 3, /? = 1 is shown in figure 9.14 for <f> p = 0 . When 
( f )p > 0°, there are no regions of monotonicity. This figure indicates (as did figure 
7.2(b) for two dimensions) that only the grid-aligned method is monotone at these 
subsonic conditions. However, it is believed that this constraint, as well as the 
constraints imposed upon supersonic flows, can be relaxed somewhat in practice 
in an effort to reduce spurious oscillations near discontinuities while still maintain- 
ing much of the sharper resolution afforded by the 5-wave model. Although an 
empirical limiting method has not been devised for three-dimensional flow due to 
its inherent complexity, it is believed that a successful method could be patterned 
much the same as the method currently employed for two dimensions. 
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Figure 9.12: Concluded 



Figure 9.13: Monotonicity for M = 3, /? - 1? Allowable n d for 




CHAPTER 10 


THREE-DIMENSIONAL EULER RESULTS 


All three-dimensional results are obtained using the implicit approximate- 
factorization algorithm CFL3D [32]. Only the flux function is varied in order to 
obtain either grid-aligned or grid-independent results. The CFL numbers chosen 
for stability for the 5-wave model are based on the theoretical analysis of Chapter 
9. No wave-propagation-direction limiting procedure has been developed for the 
three-dimensional 5- wave model for improving its monotonicity properties. 

10.1 Ramp Flow in a Channel 

The geometry for the ramp flow in a channel with inflow Mach number of 2.8 
is shown in figure 10.1. The domain is 1.8 units long by 0.72 units high by 0.5 
units wide. An oblique ramp at the lower wall starts at 0.2 units with an angle of 
12° at one side and at 0.5 units with an angle of 18° at the other side. This case 
was first computed by Parpia [33]. 

Computations are performed on a 41 x 17 x 17 grid. Two i = constant planes 
and one k = constant plane from this grid are shown in figures 10.2(a) through 
(c). Simple reflection boundary conditions are applied at the four walls, inflow is 
taken as freestream, and outflow conditions are obtained from the interior of the 
grid using uniform first-order extrapolation. 
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a) i = 1 Plane 



b) i = 17 Plane 



c) All k = Constant Planes 
Figure 10.2: 2-D Sections of 41 x 17 x 17 Grid 
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d) j = 15 Plane 

Figure 10.4: Ramp Flow, Pressure Contours, First-Order, 5- Wave 
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f) Pressure Along k = Constant Lines in i — 1 Plane 


Figure 10.5: Concluded 




d) j = 15 Plane 

Figure 10.6: Ramp Flow, Pressure Contours, Second-Order, 5- Wave 
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f) Pressure Along fc = Constant Lines m t 1 Plane 


Figure 10.6: Concluded 
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First-order computational results using the grid-aligned method are given 
in figure 10.3. Pressure contours are given in three i = constant planes, one j = 
constant plane, and one k = constant plane. Also, nondimensional pressure values 
are given along three k = constant lines in the i — 1 plane in figure 10.3(f). Results 
of the 5-wave model are obtained by restarting the grid-aligned results with the 
wave-propagation directions frozen. These 5- wave model results are given in figure 
10.4. Shock waves are resolved with far fewer interior points using this method. 

Second-order results using the grid-aligned model and the 5-wave model are 
given in figures 10.5 and 10.6, respectively. In this case, the 5-wave model does 
not improve the shock resolution to any noticeable extent over the grid-aligned 
model. Also, results are slightly more oscillatory. Hence for an oblique shock-type 
flow in three dimensions, as in two dimensions, the grid-independent model does 
not seem to be worth the additional effort when second-order spatial accuracy is 
employed. 


10.2 Oblique Supersonic Shear 


An oblique supersonic shear case is computed within a cube 1.8 units on a 
side using a Cartesian 17 x 17 x 17 mesh. For the particular case considered, the 
velocity undergoes both an increase in magnitude as well as a rotation through 
the shear layer. The shear layer itself is assumed to lie along one diagonal of the 
cube, as shown in figure 10.7. Below the layer the velocity components are u = 3, 
v = 3, and w = 3, while above the layer they are u = 4, v = 2, and w = 4. There 
is one transition cell where u = 3.5, v = 2.5, and w = 3.5. There is no pressure 
or density change across the shear layer. The boundary conditions are applied as 
follows: at inflow the exact values are specified, and at outflow the variables are 
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obtained from the interior using uniform first-order extrapolation. 

When /3 is not limited to be greater than 0.05, the 5-wave model preserves the 
exact solution when the initial condition is the exact solution. This exact solution 
is presented in figures 10.8 through 10.10 at an i = constant, a j = constant, and 
a k = constant plane, respectively. Shown are in-plane Mach number contours 
and velocity vectors (the velocity component out-of- plane is ignored). 

A first-order computation using the grid-aligned model gives the results shown 
in figures 10.11 through 10.13. The shear layer is seen to spread a significant 
amount through the domain. When restarted from the grid-aligned solution, the 
5-wave model (with f3 now limited to be no less than 0.05 to improve stability) 
gives the results shown in figures 10.14 through 10.16. For this case the wave- 
propagation directions are computed once at restart, then remain frozen for the 
remainder of the computation. The shear layer is preserved with relatively few 
interior points using the 5- wave model. Although not shown, the pressure field is 
computed in error from the exact solution (of no pressure change at all through 
the shear layer) by about 23% using the grid-aligned model, while the 5-wave 
model solution is only about 5% in error. 

In-plane Mach contours from a second-order computation using the grid- 
aligned model are given in figure 10.17, while second-order 5- wave model results, 
restarted from the first-order grid-aligned solution with the wave-propagation di- 
rections frozen, are given in figure 10.18. Even with second-order spatial accuracy 
the 5-wave model gives significantly higher shear-layer resolution than the grid- 
aligned model. The pressure field, not shown, is also computed more accurately. 
It is about 5% in error from the exact solution, while the grid-aligned result is 

about 10% in error. 
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Figure 10.7: Shear Layer Orientation Within Cube 
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b) In-Plane Velocity Vectors 


Figure 10.8: Oblique Shear, Exact Solution, * - 9 Plane 
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a) In-Plane Mach Contours b) In-Plane Velocity Vectors 


Figure 10.9: Oblique Shear, Exact Solution, j = 9 Plane 
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a) In-Plane Mach Contours b) In-Plane Velocity Vectors 


Figure 10.10: Oblique Shear, Exact Solution, fc = 9 Plane 


197 





a) In-Plane Mach Contours b) In-Plane Velocity Vectors 


Figure 10.11: Oblique Shear, First-Order, Grid- Aligned, i = 9 Plane 




a) In-Plane Mach Contours b) In-Plane Velocity Vectors 


Figure 10.12: Oblique Shear, First-Order, Grid-Aligned, j — 9 Plane 
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Oblique Shear, First-Order, Grid- Aligned, k = 9 Plane 
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Figure 10.14: Oblique Shear, First-Order, 5-Wave, i — 9 Plane 
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a) In-Plane Mach Contours b) In-Plane Velocity Vectors 

Figure 10.15: Oblique Shear, First-Order, 5- Wave, j = 9 Plane 
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a) In-Plane Mach Contours b) In-Plane Velocity Vectors 

Figure 10.16: Oblique Shear, First-Order, 5- Wave, fc = 9 Plane 


a) i = 9 Plane 


b) j = 9 Plane 


Figure 10.17: 



c) k = 9 Plane 


Oblique Shear, In-Plane Mach Contours, Second-Order, 
Grid-Aligned 
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a) i = 9 Plane 


b) j = 9 Plane 



c) A: = 9 Plane 


Figure 10.18: Oblique Shear, In-Plane Mach Contours, Second-Order, 

5- Wave 
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10.3 Transonic Flow Over a Wing 


Transonic flow is computed over the F-5 wing at the conditions M — 0.95, 
a = 0°. This wing has a root chord to tip chord to wing span ratio of 1.027 : 0.316 
: 1.0. A 97 x 17 x 17 C-H grid is used, with maximum outer boundary extent of 
about 20 wing spans, and average normal spacing at the body of 0.0118 wing spans. 
The grid extends about one wing span past the wing tip, and there are 560 cells 
on the wing surface. A cut-a-way view of this grid, with the wing delineated by a 
heavy line, is shown in figure 10.19. Standard reflection boundary conditions are 
employed at the wing surface, and subsonic inflow/outflow boundary conditions 
(as described in section 8.1.5) are employed at the farfield boundaries. 

Pressure contours at three t = constant sections over the wing and pressure 
contours on the upper surface of the wing for a first-order computation using 
the grid-aligned model are shown in figure 10.20. Pressure values are plotted in 
increments of 0.02. Since the shock waves over the wing are more or less aligned 
with grid faces, the grid-aligned method captures these very sharply. Entropy 
contours at the same locations (plotted in increments of 0.002) are shown in figure 
10 . 21 . 

The 5-wave model is run with first-order spatial accuracy by restarting the 
grid-aligned result with the wave-propagation directions frozen. Pressure contours 
are given in figure 10.22 and entropy contours are given in figure 10.23. The 5-wave 
model does not improve the shock resolution for this case. In fact, although the 
shock resolution is similar for both models over much of the wing, the 5-wave model 
actually shows worse resolution near the wing root. The reason for this is most 
likely due to the current strategy for freezing (3 and the propagation directions 
after the first iteration. Since the shocks are nearly aligned with the grid faces, 
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the grid-aligned model is already optimum for capturing them accurately. Any 
small difference in the wave-propagation direction and/or any numerical value of 
0 lower than 1.0 near the shock waves contributes to a degradation of the result. 

However, similar to the two-dimensional airfoil results, the 5-wave model also 
yields lower entropy generation over the wing than the grid-aligned model. (An 
“exact” solution would give no entropy generation except through the shocks.) 
The maximum values generated on the wing upper surface (figures 10.21(d) and 
10.23(d)) are 0.0393 for grid-aligned and 0.0284 for 5-wave. 

As this wing case demonstrates, the 5-wave model cannot give improved 
shock-resolution over the grid-aligned method when shock waves are aligned with 
grid-faces. However, as in two dimensions, it does produce lower numerically- 
generated entropy values near the body surface when first-order spatial accuracy 
is employed. Second-order computations are not performed for this case. However, 
it is expected that results would follow the trends exhibited for two-dimensional 
airfoil flow, and the 5-wave model would show no discernible advantage over the 

grid-aligned model. 
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Figure 10.19: Cut-a-Way View of F-5 Wing 97 x 17 x 17 C-H Grid 
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c) i = 9 Plane 



d) Wing Upper Surface 


Figure 10.20: Concluded 
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Figure 10.22: Concluded 




CHAPTER 11 


CONCLUSIONS 


The main object of the current research is to introduce and explore a grid- 
independent approximate Riemann solver for use with the Euler and Navier-Stokes 
equations. The primary reason for attempting a grid-independent method is to 
model the physics of the true two or three-dimensional flow more accurately. Most 
current multidimensional flow solvers that employ a Riemann solver implement 
it in a “direction-split” manner, i.e. one-dimensional theory is applied in each 
grid direction separately. In reality, however, the waves can travel in infinitely 
many directions. Constraining them to the grid directions is inconsistent with the 
physics of the flow and can result in improper interpretation of waves that are not 
aligned with the grid. 

The current grid-independent model utilizes the velocity-difference direction 
as the direction of propagation of most of the waves, with an additional shear wave 
propagating 90° to this. Use of the velocity-difference direction as the dominant 
direction allows for more accurate interpertation of both shock and shear waves 
when they lie oblique to the grid. The direction is generally frozen during a 
computation to eliminate nonlinear feedback in the solution and aid convergence. 

Simple left and right states obtained by interpolation along grid lines are used 
at each interface as the initial conditions for a Riemann problem. The difference 
in states is decomposed into a series of five waves, and the standard upwind flux 
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formula is employed to determine the flux at each face. It is shown in the derivation 
of the 5-wave model that the standard flux formula is the “best” within a certain 
class of formulas at reproducing the flux given by two rotated Riemann solvers. 

A stability analysis of the 5-wave model in conjunction with the two-dimen- 
sional Euler equations, advanced in time using explicit time-marching, shows that 
some modes of the Fourier footprint contain eigenvalues that lie on the imaginary 
axis. Hence the stability boundary of the time-marching scheme must include a 
finite portion of the imaginary axis as well. Since only multistage schemes can 
satisfy this requirement, a four-stage scheme is chosen. It is stable up to CFL 
numbers of about 1.75 for first-order spatial accuracy and 0.87 for second-order, 
although in practice these restrictions can generally be relaxed somewhat. 

An implicit approximate-factorization algorithm is shown to be stable only 
up to a CFL number of about 2.5 for first-order and 0.3 for second-order when 
the 5-wave model is employed on the right-hand side and first-order accurate 
grid-aligned approximate Jacobians are employed on the left-hand side. Again 
these CFL numbers are somewhat overrestrictive for most practical applications. 
When a grid-independent left-hand side approximate Jacobian is employed, first- 
order computations can be made unconditionally stable. Using a linearized anal- 
ysis, second-order computations are shown to be unconditionally stable if block - 
pentadiagonal systems (as opposed to tridiagonal systems) are solved with appro- 
priately chosen grid-independent approximate Jacobians. Maximum CFL num- 
bers of about 100 can be attained in practice with this strategy, although the 
optimum CFL number for convergence generally lies between 2 and 6. 

The monotonicity of the 5-wave model in conjunction with the two-dimen- 
sional Euler equations has also been investigated. Using a linearized analysis, it is 
shown that strictly only the grid-aligned first-order method is monotone. However, 
it is shown that the oscillations that occur near discontinuities can be reduced in 
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magnitude to acceptable levels for a wide variety of problems by limiting the wave- 
propagation directions of the 5-wave model according to a strategy suggested by 
the linearized analysis. This gain in monotonicity does result in the loss of some 
of the high resolution of oblique shock and shear waves, although the resolution 
is still generally greater than results using the grid-aligned scheme. 

Two-dimensional results obtained using the 5-wave model indicate that the 
method is usually worthwhile for first-order computations. Flows with oblique 
shock or shear waves are captured more sharply than results using the grid-aligned 
method, even when the 0^-limiter is employed to reduce oscillations. Also, sub- 
sonic Euler airfoil flow computations show significantly reduced entropy generation 
over the airfoil surface, resulting in better drag prediction. The 5-wave model is 
about 1.2 times more costly per iteration than the grid-aligned model when no 
^-limiting is employed, and about 1.5 times more costly when it is. Hence the 
cost penalty is not severe considering the increased accuracy of the first-order 
solutions. 

When second-order spatial accuracy is employed, the small increase in accu- 
racy attained by the 5-wave model is generally not worth the added drawbacks. In 
particular, oblique shock waves are resolved only slightly more sharply than when 
using the grid-aligned method, and there is very little decrease in the numerical 
entropy generation for a subsonic airfoil computation. These facts, taken in com- 
bination with a propensity for increased oscillatory behavior near discontinuities, 
makes the 5- wave model an unattractive alternative to the grid-aligned model in 
general. The only notable exceptions to this conclusion are found in the pure 
supersonic oblique shear wave computation and the separated viscous airfoil flow 
computation. In the first case, the 5-wave model still gives significant improved 
resolution of the shear wave over results using the grid-aligned model. In the 
second case, numerical pressure oscillations evident in the separated region over 
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the airfoil upper surface using the grid-aligned model are dramatically reduced 
by the 5-wave model, which properly interprets the oblique shear waves present 
in the flowfield. Hence in specialized cases involving oblique shear waves, the 5- 
wave model in a second-order accurate computation can be worth the additional 
expense. 

In the arena of second-order viscous flow computations, it should be noted 
that on very coarse grids the 5- wave model can do a much worse job than the grid- 
aligned model at resolving the boundary layer (for flat-plate flow) and predicting 
surface pressure values (for airfoil flow). The two models approach the same results 
as the mesh is refined, however, so as long as a fine enough grid is employed the 
5-wave model will be accurate in these areas. This problem is particularly evident 
in the viscous separated airfoil flow case, where the 5-wave model result over the 
coarsest (65 x 25) grid looks better than the grid-aligned result, since the pressure 
distortions in the flowfield have been reduced. However, the lift predicted by 
the 5-wave model is about 120% too high due to its inaccurate surface pressure 
predictions. The grid-aligned model, plagued as it is with pressure deviations in 
the separated region above the airfoil, is only high in its lift prediction by about 
8 %. 

The 5-wave model is also extended to three dimensions. The velocity-differ- 
ence direction is again chosen as the primary wave propagation direction, and the 
difference in states can still be described by a series of five waves. One of these 
waves is a shear wave which is assumed to propagate normal to the plane spanned 
by Vl and Vr. Hence the velocity of this shear wave is identically zero when the 
wave propagation directions are not frozen. Oblique shear waves through which 
the velocity undergoes a change in magnitude and/or undergoes a rotation can be 
captured sharply by this method. 

A stability analysis of the 5-wave model in three dimensions gives results 
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similar to those of the analysis in two dimensions, only with slightly lower stability 
limits. An explicit four-stage time marching scheme is stable up to CFL numbers 
of about 1.55 for first-order spatial accuracy and 0.77 for second-order. An implicit 
approximate-factorization algorithm with a first-order grid-aligned approximate 
left-hand side is stable up to a CFL number of about 1.5 for first-order and 0.3 
for second-order. 

The monotonicity properties in three dimensions are also investigated using a 
linearized analysis. Again, in a strict sense, only the grid-aligned first-order model 
turns out to be monotone. It is believed that a relaxation of this restriction on 
the wave-propagation direction similar to that in two dimensions can be employed 
in order to improve the monotonicity properties of the 5-wave model, but this 
has not been attempted in this study due to the complexity inherent in three 
dimensions. In any case, the extra degree of freedom in three dimensions seems 
to relieve some of the oscillation problems present in two dimensions when the 
5-wave model is used. 

Conclusions from three-dimensional Euler test cases run along similar lines to 
those from two dimensions. For a ramp flow with oblique shock waves, the 5-wave 
model improves the shock resolution considerably in a first-order computation, 
making the model worth the effort and additional expense. However, a second- 
order computation shows very little difference from a grid-aligned computation. 
This fact, taken in conjunction with the fact that the 5-wave model costs more and 
yields slightly increased oscillatory behavior over the grid-aligned method, makes 
the 5-wave model impractical for use with this type of problem using second-order 
spatial accuracy. 

An oblique supersonic shear wave, on the other hand, is resolved better using 
the 5-wave model for both first and second-order accurate computations. This is 
similar to the result in two dimensions. Extrapolating from the two-dimensional 
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results, then, it is likely that the 5-wave model would also do well for a three- 
dimensional separated viscous flow computation similar to the NACA 0012 airfoil 
case. Three-dimensional viscous flow computations were not attempted in this 
study. 

A first-order Euler computation performed over an F-5 wing using the 5-wave 
model produces lower entropy values over most of the the wing surface than the 
grid-aligned model. Since the shock waves for this case are aligned with grid 
faces for the most part, the 5-wave model does not improve their resolution over 
the grid-aligned model. Second-order computations on this configuration were 
not performed. However, results are expected to be similar in nature to the two- 
dimensional results, with little to no improvement in numerical entropy generation 
using the 5- wave model over the grid-aligned model. 

Overall, for both two and three-dimensional computations, the 5-wave model 
is seen to be worth the additional effort only for first-order spatially-accurate 
computations, or computations involving only oblique shear waves. In general, 
inviscid flows with shock waves and/or flows with no shock or shear waves at all 
are better computed using the grid-aligned wave model when accuracy greater 
than first-order is desired. The additional expense and oscillation-prone nature of 
the 5-wave model makes it unattractive for use in such cases. 

Despite its shortcomings, the 5-wave model can be thought of as a step in the 
right direction” toward modeling the multidimensional flow physics of the Euler 
and Navier-Stokes equations more accurately than the grid-aligned model. Future 
work may best be focused on removing the inconsistency in the model (discussed 
in section 5.4) that does not allow it to reduce to the grid-aligned model when 
Q' d = e g unless 0 = 1. Although theoretically it seems that this would be a 
desirable property to meet, it is not known what the ramifications of not meeting 
it are. Possible avenues of investigation toward eliminating the inconsistency, 
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while still maintaining sharp resolution of oblique shear waves, might include 
experimentation with different types and directions of waves as well as varying 
the flux formula itself. Further efforts toward improving the grid-independent flux 
function might focus on a different type of limiting to improve the monotonicity 
properties of the method. Perhaps by limiting the left and right input states 
to the Riemann problem, rather than the wave-propagation directions, greater 
resolution of oblique waves might be maintained for both first and second-order 
accurate computations. Also, the extension to limiting in three dimensions might 
be more straightforward than the current angle-limiting procedure. 

The current research has purposefully begun with the assumption that only 
the left and right states at each interface, obtained by interpolation along grid 
lines, are known at the beginning of the flux function computation. It is then 
up to the model to make the most of these input states, and infer from them 
the types and directions of waves likely to be present. This assumption was 
made to keep the cost of the grid-independent flux function as low as possible. 
By ignoring surrounding flowfield data, no complex interpolations are required 
to obtain gradients in non-grid-oriented directions. However, perhaps this initial 
assumption was unrealistic, and the inclusion of some multidimensional input 
states might be helpful or even necessary to improve the robustness of the model, 
particularly for second and higher-order accurate computations. 
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APPENDIX A 

DERIVATION OF (3 FOR THE MINIMUM-PATHLENGTH 
AND MINIMUM- AREA MODELS 


A.l Minimum-Pathlength Model 

As discussed in Chapter 5, the right state R lies inside the acoustic cone 
emanating from the origin when 

(Ap) 2 > [pa(Aucos6 d + Ausin^j)] 2 , (A.l) 

and the difference in states is described by two acoustic waves and an entropy 
wave. When 0 d is frozen as 6' d , then the method is modified as follows. If the 
projection (R') of the right state R onto the 0^-plane lies inside the acoustic cone, 
then two 8' d acoustic waves and an entropy wave are used, along with the necessary 
O', shear wave to connect R with R'. The condition becomes: 

a 

(A p) 2 > [paAq' d ] 2 , ( A.2) 

where A q* d is defined as 

A q' d = Aucos^ + Arsing. (A-3) 

The wavestrengths must be determined which satisfy 

4 

aw = 

jb=i 


(A.4) 
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where the +9' d acoustic, -9 d acoustic, entropy, and 9 d shear are given by: 

i t 


Pi - 
P 2 = 


CL CL , A 2 

1 , — cos# d , — sin0 d , a 
P P 


a 


1, — ^cos 9 ' d , — ; sin 9 d ,a 2 
P P 

P 3 = [i,o,o,o] T 

. T 

Pd = |0,-Tsin^,TCOS^,0 

P P 


(A.5) 


n = 


(A.6) 


These strengths turn out to be: 

& + 

# - ^ x 
£ (a 2 Ap- Ap) 

L | (— Atisinflj + Aucos^j) J 

Similarly, the condition for which one Q' d acoustic, a ( 9 d + j") shear, and an 
entropy wave (along with the necessary 9 d shear wave) are used is 


(Ap) 2 < [paAq' d ] 2 . 


(A.7) 


This indicates that the projection of the right state R onto the 0^-plane lies outside 
the acoustic cone. Again wavestrengths must be determined such that (A.4) is 
satisfied, but this time the waves are 


Pi = 


P, = 


-iT 


1, ±TCOS0^, ±^sin^,a 2 

P P 

0, -TCOsfl^-^sinfl^O 

P P 


(A.8) 


Ps = [1,0,0,0] 

p 4 = 


0, — ^sinflj, — cosflj^O 
P P 


where the + or — acoustic wave is chosen depending upon which minimizes the 
pathlength in (Au, Aw, Ap)-space. The strengths turn out to be: 


n = 


A£ 


±£? - 

^7 (a 2 Ap — Ap) 

|_ 4 (— Ausin#^ + Avcosfl^) J 


(A.9) 
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The equations (A. 6) and (A.9) can be combined as the 5-wave model described 
in Chapter 5, with strengths 

&+0&*<i'* 

& - P&Wa 
(0 - l)|A?i 
J, (a 2 Ap - Ap) 

,4 (— Ausin#^ + Avcostf^). 

for the +d' d acoustic, -9' d acoustic, (6' d + f ) shear, entropy, and 9' d shear, respec- 
tively. The variable (3 is given by 


n = 


(A.10) 


0 = min [|(| , 1] , 


(A.ll) 


where £ is defined as 

A p/(pi) ( A .12) 

A?i 

When (A.2) is true, 

ICI > 1. (A- 13 ) 

so (3 — \ and the strengths (A.6) are recovered. When (A.7) is true, (3 is assigned 
the value of C- If Ap and Aq d both have the same sign: 




a £ 


n 2 — o 

A, - 4? - ?A,i 


(A. 14) 


and (A.9) with the + acoustic wave is recovered. If Ap and Aq d have opposite 
signs: 


= 0 

si, = 4? 


a* 


(A. 15) 


fla — 


Ap p 

' A2 


-^Aq' d 


a* a 

and (A.9) with the —acoustic wave is recovered. In both cases and 0$ are the 
same, as given by (A.10). With this method, the proper acoustic wave is always 
chosen such that the (8 d + f ) shear wavestrength fi 3 is minimized. 
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Notice that the wave patterns change in a continuous fashion as £ goes 
through 1.0, as shown in figure A.l. 


A. 2 Minimum-Area Model 


As discussed in Chapter 5, the minimum-area model picks the path in (Au, 
Au, Ap)-space that lies, in a sense, closest to the straight line connecting L and 
R. It does this by minimizing the area between the paths and the straight line. 
When 9 d is frozen as 9 f d , this criterion is modified slightly for simplicity. A 9' d 
shear wave is used to connect the right state R with R\ its projection onto the 
0^-plane. Then +9' d acoustic, -0' d acoustic, and (9 f d + f ) shear waves are chosen 
which minimize the area between the wave paths and the straight line connecting 
L and R' in ( Au, At?, Ap)-space. (The entropy wave, as always, is present but 
only contributes to changes in density. Its strength is always Q 5 = Ap - A p/a 2 .) 

When inequality (A. 2) is true, two acoustics and no (9 d + f ) shear wave 
minimize the area. This is identical to the minimum-pathlength model and hence 
corresponds to /3 = 1 in (A. 10). When (A. 7) is true, some combination of all 
three waves minimizes the area. The order in which the waves are taken makes 
a difference in the outcome. In the present derivation they are always taken 
such that the first wave (from L) is an acoustic that points away from R in 
(Au, Au, Ap)-space. The second wave is the ( 9 d + y) shear, and the third wave 
is the remaining acoustic. This strategy is illustrated in figure A. 2. The type 
of acoustic that leaves L depends on the sign of A p and the sign of A q d . If A p 
and A q d have the same sign, a —9 d acoustic leaves L, while if Ap and A q d have 
opposite signs, a -\-9 d acoustic leaves L. 

In this derivation it is assumed that the waves have strengths given by (A. 10), 



223 



224 


where /? is unknown. An expression for /3 is then determined that gives the 
minimum area. Figure A. 3 is the same as figure A. 2 with some relevant distances 
labeled. It is assumed for now that A p and A^ are both positive. Hence the first 
wave to leave L is a — 0 d acoustic. It turns out that the end result for [3 is the 
same regardless of the signs of Ap and A q' d . 

From the geometry of figure A.3, the areas of the two triangles in the figure 


are 


Vi = -^( Axj + Ax s )Ap! 

z 

V 2 = ^Ax 4 Ap 2 . 


(A.16) 


Since it is known (from (A.9)) that an acoustic wave which causes a change in 
pressure Ap has an associated strength Ap/ a 2 , Api can be found since the strength 
of the —acoustic wave is known: 


Cl 


— ac 


Ap 

2fi2 




(A.17) 


Therefore 

Api = ~ - ^0paAq' d . (A. 18) 

The changes in primitive variables caused by a ±0^ acoustic wave are: 

iL = ± Su = ± t-. (A. 19) 

pa 2 acost^ asm0' d 

From geometry, this can also be written as 

^ = ±(£ucos0^ + £usin<%). (A. 20) 

pa 

This means that the distance Axi is equal to Api/(pa), or 

Axi = ~ ~/3A q d . ( A - 21 ) 

2 pa 2 

Furthermore, 

Ap 2 = Ap - Api = -y + ^(3paAq d 


(A.22) 


225 


A x ; = ^ = £? + 

pa Zpa Z 

The distance Ax 3 is found using the riile of similar triangles 

Ax 3 _ Api 
A q' d A p ’ 


or 


Ax 3 — — Ag[j ~~ 2A * 


(A. 23) 


(A.24) 


(A. 25) 


Then Ax 4 is simply 

Ax 4 = Ag^ - Ax 2 - Ax 3 


1 A 1 A ? , 
= 2 A, ‘ i “2M + 




(A. 26) 


Now, the total area V T = V 1 + V 2 is found by plugging in the appropriate expres- 
sions into (A. 16). The final result is 

V T = iA P A,i - p {iApA,i} +0> {^(A ? i)’} • (A.27) 

To find the value of (3 for which V T is minimized, solve for /? when dV T /d(3 = 0: 

P _ f *m y . (A. 28) 

Also, d 2 V T /dp 2 = p 2 a 2 (Aq' d y/(2Ap) > 0 for Ap > 0 and A q' d > 0, so this (3 
indeed represents the minimum area. 

It turns out that the expression for total area is given by (A.27) when Ap 
and A q' d have the same sign, but is given by the negative of (A.27) when they 
have opposite signs. Therefore (A.28) is the correct expression for R' outside of 
the acoustic cone regardless of the sign of Ap and Aq d . It gives strengths of 

n 1 = ^?+ V 


Ci ■> — 


2a 2 2pa 3 Aq , d 
Ap Ap 2 
2 a 2 2pa 3 Aq d 

ft > = 4£r - f a,; 

pa 3 Aq d a 


(A.29) 
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for the +0' d acoustic, -d' d acoustic, and (0' d + f ) shear, respectively. 

The results for /? for the two cases of R' outside or inside the acoustic cone 
can be combined as: 

/? = min [< 2 ,l] . ( A.30) 

where again ( is given by (A. 12). The wave patterns represented by this model 
change in a continuous fashion as ( goes through 1.0, as shown in figure A. 4. 
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APPENDIX B 

DERIVATION OF THE COMBINED SHEAR- WAVE VECTOR 
FOR THE 3-D EULER EQUATIONS 


The eigenvalues of the linearized system of three-dimensional Euler equations 
can be written as 

Ai = Qg + 0. 

X 2 = q g — a (B-l) 

A3 = A4 = 

where q g is defined by (9.28). The eigenvector for the conserved-variable form of 
the equations associated with Ai is given by Ri in (9.23), while the eigenvector 
associated with A 2 is given by R 2 in (9.24). 

There are three linearly independent eigenvectors associated with the eigen- 
value q g that form a basis for its eigenspace. One, an entropy wave, is denoted by 
equation (9.26). The other two are shear waves that cause a change in the velocity 
with no change in density or pressure. Any two shear vectors that are orthogonal 
when written in primitive- variable form and, cause a change in velocity normal to 
the fig -direction can be chosen. It turns out that the end result for the combined 
shear vectors is the same regardless of the choice. One set is 


R 


a hear 1 


0 

— asin0 g 
a cosQ g 
0 

_a(— usin0 fl +t)cos0 9 ). 


(B.2) 
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R-shear2 


0 

—acos9 g siml)g 

—asindgsintpg , 

a cosi/} g 

.a(—ucos9 g s\mpg — vsm9 g simj> g + iocosi /> 9 ). 


(B.3) 


where the angles 9 g and define the grid-face normal direction as pictured in 


figure 9.1. 

The strengths of each of the five waves can be obtained by solving 


AU = £fi fe R* (B.4) 

k 

for the f l k . Only the shear-wave strengths are necessary in the present derivation. 
They are given by: 

fishearl = T(-Ausin0 fl + AuCOsflj) ' 

o (B.5) 

ft shear2 = t (- A ttcos0 9 sin0 ff - AusinfljSinV’,, + Awcos if} g ). 
a 

The combined action of the two shear waves in the linearized system can be 
replaced by a single shear wave whose strength and direction are a function of the 
difference in states. An expression for this new shear vector can be found from: 


a 


newUnew — OahearlR'shearl T flshear2 R'shear 2 ■ 


(B.6) 


Combining terms, one obtains: 


OnewR-new — 


~p{ An), 

— p(Ar 2 ) 9 
p(Ar 3 ) 9 

L-pfufAr^g + u(Ar 2 ) 9 to(Ar 3 ) 9 } J 


where the (Arj) s ’s are defined by 

(Ari), = ((c.)J - 1) An + (c x )g{c y ) g Av + (c z ) g (c z ) g Aw 
(Ar 2 ) 9 = (c s ) 9 (c x ) s Au + ((c v )g - l)Aw + (c y )g(c z )gAw 
(Ar 3 ) 9 = [c z ) g (c x ) g Au 4- (c z ) 9 (c y ) 9 Av + ((c z ) 9 — l) A uj. 


(B.7) 


(B.8) 
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The c x , c y , and c z are defined in spherical coordinates in (9.20). In order that the 
vector R new be somewhat similar in form to the two vectors from which it was 
derived, fi new is taken to be p/a\ the shear vector R ne w then becomes: 


R 


new 


o 

— a( At'i ) 9 

-a(Ar 2 ) 9 
— a(Ar 3 ) 9 

. — d{u(Ari) 9 + v(Ar 2 ) 9 + d)(Ar 3 ) 9 }. 


(B.9) 


Notice that A q g = (c x ) g Au + (c y ) g Av + (c z ) g Aw represents the component 
of the velocity-difference AV^ in the wave-propagation direction (normal to the 
grid face), so (B.8) can be rewritten as 

(^ r i)a = Aq g (c x ) g — Au 

(A r 2 ) g = A q g (c y ) g - Av (B.10) 

(A r 3 )g = Aq g (c y )g — Aw. 

Hence each of these terms represents the difference between the x, y, or z com- 
ponent of the velocity difference along the grid-normal and the corresponding 
component of the whole velocity difference. 
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